### File    : corstange-jop-syria-replication.r
##  Author  : Daniel Corstange (daniel.corstange@columbia.edu)
##  Created : Wed 23 Jun 2021 --- 03:42 PM
##  Notes   : data analysis for Assad/jihadi list experiment

### read/process data
#### preliminaries

## load various libraries
library(haven)
require(MASS)
library(texreg)
library(ggplot2)
library(dplyr)
library(tidyr)
library(lubridate)
library(stringr)
library(xtable)
library(egg)
library(list)

#### raw data

Syr17 <- read_spss("Syria-2017.sav")
Syr15 <- read_spss("Syria-2015.sav")

#### process data
##### background covariates
###### 2017

BG17 <-
    transmute(
        Syr17,
        ## metadata
        rid    = 170000 + 1:nrow(Syr17),
        year   = 2017,
        ## months as refugee
        ## Q14b = year, Q14a = month
        ymd    = paste(Q14b, Q14a, "1", sep = "-"),
        months =
            as.numeric(
                as.duration(interval(ymd, "2017-07-01")),
                "months"
            ),
        province_syr =
            factor(as_factor(Q13)),
        province_leb =
            factor(as_factor(Mohafaza)),
        district =
            recode(factor(as_factor(Caza)),
                   "Miniyé-Danniyé" = "Minieh-Dinnieh",
                   "Nabatiyé" = "Nabatieh",
                   "West Békaa" = "West Bekaa",
                   "Tripoli (Trablous)" = "Tripoli",
                   "El Metn" = "Metn",
                   "Aakkar" = "Akkar"
                   ),
        ## basic demographics
        female = as.logical(as.numeric(Q5)),
        male   = !female,
        age    = Q6,
        age30p = ifelse(age >= 30, TRUE, FALSE),
        young  = !age30p,
        ## room density
        syr_room_density = Q8 / (Q9 + 1),
        leb_room_density = Q10 / (Q11 + 1),
        syr_fam_size = Q8,
        leb_fam_size = Q10,
        srd2p = ifelse(syr_room_density >= 2, TRUE, FALSE),
        lrd2p = ifelse(leb_room_density >= 2, TRUE, FALSE),
        ## finished high school
        ## watch out for the backquote (`) and the leading space!
        educ_literate =
            ifelse(as_factor(Q7) != "Illiterate",
                   TRUE, FALSE),
        educ_hs = ifelse(as_factor(Q7)
                         %in% c("Finished secondary school",
                                "Some university",
                                "Bachelor`s degree",
                                "Master`s degree or higher",
                                " Vocational/Technical training"),
                         TRUE, FALSE),
        educ_college =
            ifelse(as_factor(Q7)
                   %in% c("Bachelor`s degree",
                          "Master`s degree or higher",
                          " Vocational/Technical training"),
                   TRUE, FALSE),
        ## knowledge
        know = ((as_factor(Q23a) == "Walid al-Mouallem") +
                (as_factor(Q23b) == "7 years") +
                (as_factor(Q23c) == "250") +
                (as_factor(Q23d) == "Tunisia")),
        know_bi = ifelse(know > 2, TRUE, FALSE),
        ## interest 4-point, 0-3 low-high
        interest = ifelse(Q21 <= 0, NA, abs(Q21 - 4)),
        ## interest: not interested
        interest_low = ifelse(Q21 <= 0, NA,
                       ifelse(as_factor(Q21) == "Not interested", TRUE,
                              FALSE)),
        ## somewhat/very interested in politics
        interest_bi = ifelse(as_factor(Q21)
                             %in% c("Very interested",
                                    "Somewhat interested"), TRUE,
                      ifelse(as_factor(Q21)
                             %in% c("A little interested",
                                    "Not interested"), FALSE,
                             NA)),
        ## understand politics (middle category counted as "understand")
        understand_bi =
            ifelse(Q22 <= -98, NA,
            ifelse(Q22 >= 3, TRUE,
                   FALSE)),
        ## index: add number of true items
        engaged = know_bi + interest_bi + understand_bi,
        engaged_bi =                     # dichotomize
            ifelse(engaged > 1,
                   TRUE, FALSE),
        ## community
        ## watch out for backquote (`)!
        sunni =
            ifelse(as_factor(Q16) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q16) == "Sunni", TRUE,
                   FALSE)),
        kurdish_speak =
            ifelse(as_factor(Q15b) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
            ifelse(as_factor(Q15b) == "Comfortably", TRUE,
                   FALSE)),
        minority = ifelse(!sunni | kurdish_speak, TRUE, FALSE),
        ## religiosity
        pray_daily =
            ifelse(as_factor(Q17) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q17) %in% c("Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        pray_multiple =
            ifelse(as_factor(Q17) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q17) %in% c("Several times a day"), TRUE,
                   FALSE)),
        quran_daily =
            ifelse(as_factor(Q18) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q18) %in% c("Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        quran_weekly =
            ifelse(as_factor(Q18) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q18) %in% c("Once a week",
                                         "Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        ## religion and politics
        rel_econ =
            ifelse(as_factor(Q19a) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
                   Q19a),
        rel_pol =
            ifelse(as_factor(Q19b) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
                   Q19b),
        rel_crime =
            ifelse(as_factor(Q19c) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
                   Q19c),
        rel_family =
            ifelse(as_factor(Q19d) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
                   Q19d),
        ## index: add number of items at least "somewhat important"
        relpol =                        # additive index
            ifelse(rel_econ   > 2, 1, 0) +
            ifelse(rel_pol    > 2, 1, 0) +
            ifelse(rel_crime  > 2, 1, 0) +
            ifelse(rel_family > 2, 1, 0),
        relpol_bi =                     # dichotomize
            ifelse(relpol > 2,
                   TRUE, FALSE),
        ## factional support
        fac_none =            # no faction mentioned
            ifelse(Q24a == -99 &
                   Q24b == -99 &
                   Q24c == -99 &
                   Q24d == -99 &
                   Q24e == -99 &
                   Q24f == -99,
                   TRUE, FALSE),
        fac_govt =            # if govt 1st or 2nd (most Kurds say govt 2nd)
            ifelse(Q24b %in% 1:2,
                   TRUE, FALSE),
        fac_fsa =             # fsa 1st, NOT mention Islamists
            ifelse(Q24a == 1 &
                   Q24c <  0 &
                   Q24d <  0,
                   TRUE, FALSE),
        fac_islamist =        # any Islamists mentioned
            ifelse(Q24c > 0 |
                   Q24d > 0,
                   TRUE, FALSE),
        fac_islamistd =       # domestic mentioned, foreign NOT mentioned
            ifelse(Q24c > 0 &
                   Q24d < 0,
                   TRUE, FALSE),
        fac_islamistf =       # foreign mentioned at all
            ifelse(Q24d > 0,
                   TRUE, FALSE),
        fac_oppn =
            ifelse(fac_govt == 0 &
                   fac_none == 0,
                   TRUE, FALSE),
        fac_all =
            factor(ifelse(fac_none, 1,
                   ifelse(fac_islamistf, 2,
                   ifelse(fac_islamistd, 3,
                   ifelse(fac_fsa, 4,
                   ifelse(fac_govt, 5,
                          6))))),
                   labels = c("No Preference",
                              "Foreign Islamists",
                              "Domestic Islamists",
                              "Nationalists",
                              "Government",
                              "Other")),
        ## convenience faction variables
        fac_4cat =                      # none, other, govt, oppn
            factor(ifelse(fac_all == "No Preference", 1,
                   ifelse(fac_all == "Other", 2,
                   ifelse(fac_all == "Government", 3,
                          4))),
                   labels = c("none",
                              "other",
                              "govt",
                              "oppn")),
        fac_5cat =                      # none, other, govt, fsa, islamists
            factor(ifelse(fac_all == "No Preference", 1,
                   ifelse(fac_all == "Other", 2,
                   ifelse(fac_all == "Government", 3,
                   ifelse(fac_all == "Nationalists", 4,
                          5)))),
                   labels = c("none",
                              "other",
                              "govt",
                              "fsa",
                              "islamists"))
    )

###### 2015

BG15 <-
    transmute(
        Syr15,
        ## metadata
        rid    = 150000 + 1:nrow(Syr15),
        year   = 2015,
        province_syr =
            factor(as_factor(Q15)),
        province_leb =
            factor(as_factor(MOHFZ)),
        district =
            recode(factor(as_factor(CAZA)),
                   "Minieh & Dinnieh" = "Minieh-Dinnieh",
                   "Tyre" = "Sour"
                   ),
        ## basic demographics
        female = as.logical(as.numeric(Q5)),
        male   = !female,
        age    = Q7,
        age30p = ifelse(age >= 30, TRUE, FALSE),
        young  = !age30p,
        ## months as refugee
        ## Q16b = year, Q16a = month
        ymd    = paste(Q16b, Q16a, "1", sep = "-"),
        months =
            as.numeric(
                as.duration(interval(ymd, "2015-05-01")),
                "months"
            ),
        ## room density
        syr_room_density = Q10 / (Q11 + 1),
        leb_room_density = Q12 / (Q13 + 1),
        syr_fam_size = Q10,
        leb_fam_size = Q12,
        srd2p = ifelse(syr_room_density >= 2, TRUE, FALSE),
        lrd2p = ifelse(leb_room_density >= 2, TRUE, FALSE),
        ## finished high school
        ## watch out for the backquote (`) and the leading space!
        educ_literate =
            ifelse(as_factor(Q8) != "Illiterate",
                   TRUE, FALSE),
        educ_hs = ifelse(as_factor(Q8)
                         %in% c("Finished secondary school",
                                "Some university",
                                "Bachelor`s degree",
                                "Master`s degree or higher",
                                "Vocational/Technical training"),
                         TRUE, FALSE),
        educ_college =
            ifelse(as_factor(Q8)
                   %in% c("Bachelor`s degree",
                          "Master`s degree or higher",
                          "Vocational/Technical training"),
                   TRUE, FALSE),
        ## knowledge
        know = ((as_factor(Q27a) == "Walid al-Mouallem") +
                (as_factor(Q27c) == "7 Years") +
                (as_factor(Q27d) == "250") +
                (as_factor(Q27e) == "Tunisia")),
        know_bi = ifelse(know > 2, TRUE, FALSE),
        ## interest 4-point, 0-3 low-high
        interest = ifelse(Q25 <= 0, NA, abs(Q25 - 4)),
        ## interest: not interested
        interest_low = ifelse(Q25 <= 0, NA,
                       ifelse(as_factor(Q25) == "Not interested", TRUE,
                              FALSE)),
        ## somewhat/very interested in politics
        interest_bi = ifelse(as_factor(Q25)
                             %in% c("Very interested",
                                    "Somewhat interested"), TRUE,
                      ifelse(as_factor(Q25)
                             %in% c("A little interested",
                                    "Not interested"), FALSE,
                             NA)),
        ## understand politics (middle category counted as "understand")
        understand_bi = ifelse(Q26 <= -98, NA,
                        ifelse(Q26 >= 3, TRUE,
                               FALSE)),
        ## index: add number of true items
        engaged = know_bi + interest_bi + understand_bi,
        engaged_bi =                     # dichotomize
            ifelse(engaged > 1,
                   TRUE, FALSE),
        ## community
        ## watch out for backquote (`)!
        sunni =
            ifelse(as_factor(Q21) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q21) == "Sunni", TRUE,
                   FALSE)),
        kurdish_speak =
            ifelse(as_factor(Q19b) %in% c("No Answer",
                                          "Don`t Know",
                                          "Not Applicable"), NA,
            ifelse(as_factor(Q19b) == "Comfortably", TRUE,
                   FALSE)),
        minority = ifelse(!sunni | kurdish_speak, TRUE, FALSE),
        ## religiosity
        pray_daily =
            ifelse(as_factor(Q22) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q22) %in% c("Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        pray_multiple =
            ifelse(as_factor(Q22) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q22) %in% c("Several times a day"), TRUE,
                   FALSE)),
        quran_daily =
            ifelse(as_factor(Q23) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q23) %in% c("Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        quran_weekly =
            ifelse(as_factor(Q23) %in% c("No Answer",
                                         "Don`t Know",
                                         "Not Applicable"), NA,
            ifelse(as_factor(Q23) %in% c("Once a week",
                                         "Once or twice a day",
                                         "Several times a day"), TRUE,
                   FALSE)),
        ## religion and politics
        rel_econ =
            ifelse(as_factor(Q24a) %in% c("No Answer",
                                          "Don`t Know",
                                          "Don`t know",
                                          "Not Applicable"), NA,
                   Q24a),
        rel_pol =
            ifelse(as_factor(Q24b) %in% c("No Answer",
                                          "Don`t Know",
                                          "Don`t know",
                                          "Not Applicable"), NA,
                   Q24b),
        rel_crime =
            ifelse(as_factor(Q24c) %in% c("No Answer",
                                          "Don`t Know",
                                          "Don`t know",
                                          "Not Applicable"), NA,
                   Q24c),
        rel_family =
            ifelse(as_factor(Q24d) %in% c("No Answer",
                                          "Don`t Know",
                                          "Don`t know",
                                          "Not Applicable"), NA,
                   Q24d),
        ## index: add number of items at least "somewhat important"
        relpol =                        # additive index
            ifelse(rel_econ   > 2, 1, 0) +
            ifelse(rel_pol    > 2, 1, 0) +
            ifelse(rel_crime  > 2, 1, 0) +
            ifelse(rel_family > 2, 1, 0),
        relpol_bi =                     # dichotomize
            ifelse(relpol > 2,
                   TRUE, FALSE),
        ## factional support
        fac_none =                # no faction mentioned
            ifelse(is.na(Q28a) &
                   is.na(Q28b) &
                   is.na(Q28c) &
                   is.na(Q28d) &
                   is.na(Q28e) &
                   is.na(Q28f),
                   TRUE, FALSE),
        fac_govt =                # if govt 1st or 2nd (most Kurds say govt 2nd)
            ifelse(is.na(Q28b) | Q28b == 3,
                   FALSE, TRUE),
        fac_fsa =                 # fsa 1st, NOT mention Islamists
            ifelse(is.na(Q28a)  |
                   Q28a != 1    |
                   !is.na(Q28c) |
                   !is.na(Q28d),
                   FALSE, TRUE),
        fac_islamist =            # any Islamists mentioned
            ifelse(is.na(Q28c) &
                   is.na(Q28d),
                   FALSE, TRUE),
        fac_islamistd =           # domestic mentioned, foreign NOT mentioned
            ifelse(!is.na(Q28c) &
                   is.na(Q28d),
                   TRUE, FALSE),
        fac_islamistf =           # foreign mentioned at all
            ifelse(is.na(Q28d),
                   FALSE, TRUE),
        fac_oppn =
            ifelse(fac_govt &
                   fac_none,
                   TRUE, FALSE),
        fac_all =
            factor(ifelse(fac_none, 1,
                   ifelse(fac_islamistf, 2,
                   ifelse(fac_islamistd, 3,
                   ifelse(fac_fsa, 4,
                   ifelse(fac_govt, 5,
                          6))))),
                   labels = c("No Preference",
                              "Foreign Islamists",
                              "Domestic Islamists",
                              "Nationalists",
                              "Government",
                              "Other")),
        ## convenience faction variables
        fac_4cat =                      # none, other, govt, oppn
            factor(ifelse(fac_all == "No Preference", 1,
                   ifelse(fac_all == "Other", 2,
                   ifelse(fac_all == "Government", 3,
                          4))),
                   labels = c("none",
                              "other",
                              "govt",
                              "oppn")),
        fac_5cat =                      # none, other, govt, fsa, islamists
            factor(ifelse(fac_all == "No Preference", 1,
                   ifelse(fac_all == "Other", 2,
                   ifelse(fac_all == "Government", 3,
                   ifelse(fac_all == "Nationalists", 4,
                          5)))),
                   labels = c("none",
                              "other",
                              "govt",
                              "fsa",
                              "islamists"))
    )

##### Assad/caliphate list
###### 2017

AC_list17 <- 
    ## keep only useful variables...
    select(
        ## ...but there are going to be some helper variables
        transmute(
            Syr17,
            ## metadata
            rid    = 170000 + 1:nrow(Syr17),
            year   = 2017,
            ## treatment groups
            ac_group =
                factor(x = Q35,
                       labels = c("g1:control-assad-good",
                                  "g2:control-caliphate-good",
                                  "g3:treat-assad-good",
                                  "g4:treat-caliphate-good")),
            ac_group_control =
                grepl("control", ac_group),
            ac_group_treat =
                grepl("treat", ac_group),
            ac_group_assad =
                grepl("assad", ac_group),
            ac_group_caliphate =
                grepl("caliphate", ac_group),
            ## convenience variables (simplifies working with tidyr)
            ac_treat_type =             # control or treatment group?
                factor(ifelse(grepl("control",
                                    ac_group),
                              1, 2),
                       labels = c("ctrl",
                                  "treat")),
            ac_assad_type =             # Assad good, or just okay?
                factor(ifelse(grepl("assad",
                                    ac_group),
                              1, 2),
                       labels = c("good",
                                  "okay")),
            ac_caliphate_type =         # caliphate good, or just okay?
                factor(ifelse(grepl("caliphate",
                                    ac_group),
                              1, 2),
                       labels = c("good",
                                  "okay")),
            ## helper variables for group 1
            ## (DK/NR is -99; not in group is -98)
            ac_parl_g1 =
                ifelse(is.na(Q36a), -98, Q36a),
            ac_no_parties_g1 =
                ifelse(is.na(Q36b), -98, Q36b),
            ac_assad_g1 =
                ifelse(is.na(Q36c), -98, Q36c),
            ac_mil_g1 =
                ifelse(is.na(Q36d), -98, Q36d),
            ac_pres_g1 =
                ifelse(is.na(Q37a), -98, Q37a),
            ac_multiparty_g1 =
                ifelse(is.na(Q37b), -98, Q37b),
            ac_caliphate_g1 =
                ifelse(is.na(Q37c), -98, Q37c),
            ac_mandate_g1 =
                ifelse(is.na(Q37d), -98, Q37d),
            ## helper variables for group 2
            ## (DK/NR is -99; not in group is -98)
            ac_parl_g2 =
                ifelse(is.na(Q39a), -98, Q39a),
            ac_no_parties_g2 =
                ifelse(is.na(Q39b), -98, Q39b),
            ac_assad_g2 =
                ifelse(is.na(Q39c), -98, Q39c),
            ac_mil_g2 =
                ifelse(is.na(Q39d), -98, Q39d),
            ac_pres_g2 =
                ifelse(is.na(Q38a), -98, Q38a),
            ac_multiparty_g2 =
                ifelse(is.na(Q38b), -98, Q38b),
            ac_caliphate_g2 =
                ifelse(is.na(Q38c), -98, Q38c),
            ac_mandate_g2 =
                ifelse(is.na(Q38d), -98, Q38d),
            ## helper variables for group 3
            ## (DK/NR is -99; not in group is -98)
            ac_assad_g3 =
                ifelse(is.na(Q40), -98, Q40),
            ac_caliphate_g3 =
                ifelse(is.na(Q41), -98, Q41),
            ## helper variables for group 4
            ac_assad_g4 =
                ifelse(is.na(Q43), -98, Q43),
            ac_caliphate_g4 =
                ifelse(is.na(Q42), -98, Q42),
            ## direct responses to each list item (groups 1 and 2)
            ac_parl_direct =
                ifelse(ac_parl_g1 == 1 |
                       ac_parl_g2 == 1, TRUE,
                ifelse(ac_parl_g1 == 0 |
                       ac_parl_g2 == 0, FALSE,
                       NA)),
            ac_no_parties_direct =
                ifelse(ac_no_parties_g1 == 1 |
                       ac_no_parties_g2 == 1, TRUE,
                ifelse(ac_no_parties_g1 == 0 |
                       ac_no_parties_g2 == 0, FALSE,
                       NA)),
            ac_assad_direct =
                ifelse(ac_assad_g1 == 1 |
                       ac_assad_g2 == 1, TRUE,
                ifelse(ac_assad_g1 == 0 |
                       ac_assad_g2 == 0, FALSE,
                       NA)),
            ac_mil_direct =
                ifelse(ac_mil_g1 == 1 |
                       ac_mil_g2 == 1, TRUE,
                ifelse(ac_mil_g1 == 0 |
                       ac_mil_g2 == 0, FALSE,
                       NA)),
            ac_pres_direct =
                ifelse(ac_pres_g1 == 1 |
                       ac_pres_g2 == 1, TRUE,
                ifelse(ac_pres_g1 == 0 |
                       ac_pres_g2 == 0, FALSE,
                       NA)),
            ac_multiparty_direct =
                ifelse(ac_multiparty_g1 == 1 |
                       ac_multiparty_g2 == 1, TRUE,
                ifelse(ac_multiparty_g1 == 0 |
                       ac_multiparty_g2 == 0, FALSE,
                       NA)),
            ac_caliphate_direct =
                ifelse(ac_caliphate_g1 == 1 |
                       ac_caliphate_g2 == 1, TRUE,
                ifelse(ac_caliphate_g1 == 0 |
                       ac_caliphate_g2 == 0, FALSE,
                       NA)),
            ac_mandate_direct =
                ifelse(ac_mandate_g1 == 1 |
                       ac_mandate_g2 == 1, TRUE,
                ifelse(ac_mandate_g1 == 0 |
                       ac_mandate_g2 == 0, FALSE,
                       NA)),
            ## list responses to the list (groups 3 and 4; add up g1 and g2)
            ac_assad_list =
                ifelse(ac_assad_g3 >= 0, ac_assad_g3,
                ifelse(ac_assad_g4 >= 0, ac_assad_g4,
                       ac_parl_direct +
                       ac_no_parties_direct +
                       ac_mil_direct)),
            ac_caliphate_list =
                ifelse(ac_caliphate_g3 >= 0, ac_caliphate_g3,
                ifelse(ac_caliphate_g4 >= 0, ac_caliphate_g4,
                       ac_pres_direct +
                       ac_multiparty_direct +
                       ac_mandate_direct))
        ),
        ## keep only the nonnhelper variables
        rid,
        year,
        ac_group,
        ac_group_control,
        ac_group_treat,
        ac_group_assad,
        ac_group_caliphate,
        ac_treat_type,
        ac_assad_type,
        ac_caliphate_type,
        ac_parl_direct,
        ac_no_parties_direct,
        ac_assad_direct,
        ac_mil_direct,
        ac_pres_direct,
        ac_multiparty_direct,
        ac_caliphate_direct,
        ac_mandate_direct,
        ac_assad_list,
        ac_caliphate_list
    )

###### 2015

AC_list15 <- 
    ## keep only useful variables...
    select(
        ## ...but there are going to be some helper variables
        transmute(
            Syr15,
            ## metadata
            rid    = 150000 + 1:nrow(Syr15),
            year   = 2015,
            ## treatment groups
            ac_group = factor(x = Q43,
                              labels = c("g1:control-assad-good",
                                         "g2:control-caliphate-good",
                                         "g3:treat-assad-good",
                                         "g4:treat-caliphate-good")),
            ac_group_control =
                grepl("control", ac_group),
            ac_group_treat =
                grepl("treat", ac_group),
            ac_group_assad =
                grepl("assad", ac_group),
            ac_group_caliphate =
                grepl("caliphate", ac_group),
            ## convenience variables (simplifies working with tidyr)
            ac_treat_type =             # control or treatment group?
                factor(ifelse(grepl("control",
                                    ac_group),
                              1, 2),
                       labels = c("ctrl",
                                  "treat")),
            ac_assad_type =             # Assad good, or just okay?
                factor(ifelse(grepl("assad",
                                    ac_group),
                              1, 2),
                       labels = c("good",
                                  "okay")),
            ac_caliphate_type =         # caliphate good, or just okay?
                factor(ifelse(grepl("caliphate",
                                    ac_group),
                              1, 2),
                       labels = c("good",
                                  "okay")),
            ## helper variables for group 1
            ## (DK/NR is -99; not in group is -98)
            ac_parl_g1 =
                ifelse(is.na(Q44a), -98, Q44a),
            ac_no_parties_g1 =
                ifelse(is.na(Q44b), -98, Q44b),
            ac_assad_g1 =
                ifelse(is.na(Q44c), -98, Q44c),
            ac_mil_g1 =
                ifelse(is.na(Q44d), -98, Q44d),
            ac_pres_g1 =
                ifelse(is.na(Q45a), -98, Q45a),
            ac_multiparty_g1 =
                ifelse(is.na(Q45b), -98, Q45b),
            ac_caliphate_g1 =
                ifelse(is.na(Q45c), -98, Q45c),
            ac_mandate_g1 =
                ifelse(is.na(Q45d), -98, Q45d),
            ## helper variables for group 2
            ## (DK/NR is -99; not in group is -98)
            ac_parl_g2 =
                ifelse(is.na(Q47a), -98, Q47a),
            ac_no_parties_g2 =
                ifelse(is.na(Q47b), -98, Q47b),
            ac_assad_g2 =
                ifelse(is.na(Q47c), -98, Q47c),
            ac_mil_g2 =
                ifelse(is.na(Q47d), -98, Q47d),
            ac_pres_g2 =
                ifelse(is.na(Q46a), -98, Q46a),
            ac_multiparty_g2 =
                ifelse(is.na(Q46b), -98, Q46b),
            ac_caliphate_g2 =
                ifelse(is.na(Q46c), -98, Q46c),
            ac_mandate_g2 =
                ifelse(is.na(Q46d), -98, Q46d),
            ## helper variables for group 3
            ## (DK/NR is -99; not in group is -98)
            ac_assad_g3 =
                ifelse(is.na(Q48z), -98, Q48z),
            ac_caliphate_g3 =
                ifelse(is.na(Q49z), -98, Q49z),
            ## helper variables for group 4
            ac_assad_g4 =
                ifelse(is.na(Q51z), -98, Q51z),
            ac_caliphate_g4 =
                ifelse(is.na(Q50z), -98, Q50z),
            ## direct responses to each list item (groups 1 and 2)
            ac_parl_direct =
                ifelse(ac_parl_g1 == 1 |
                       ac_parl_g2 == 1, TRUE,
                ifelse(ac_parl_g1 == 0 |
                       ac_parl_g2 == 0, FALSE,
                       NA)),
            ac_no_parties_direct =
                ifelse(ac_no_parties_g1 == 1 |
                       ac_no_parties_g2 == 1, TRUE,
                ifelse(ac_no_parties_g1 == 0 |
                       ac_no_parties_g2 == 0, FALSE,
                       NA)),
            ac_assad_direct =
                ifelse(ac_assad_g1 == 1 |
                       ac_assad_g2 == 1, TRUE,
                ifelse(ac_assad_g1 == 0 |
                       ac_assad_g2 == 0, FALSE,
                       NA)),
            ac_mil_direct =
                ifelse(ac_mil_g1 == 1 |
                       ac_mil_g2 == 1, TRUE,
                ifelse(ac_mil_g1 == 0 |
                       ac_mil_g2 == 0, FALSE,
                       NA)),
            ac_pres_direct =
                ifelse(ac_pres_g1 == 1 |
                       ac_pres_g2 == 1, TRUE,
                ifelse(ac_pres_g1 == 0 |
                       ac_pres_g2 == 0, FALSE,
                       NA)),
            ac_multiparty_direct =
                ifelse(ac_multiparty_g1 == 1 |
                       ac_multiparty_g2 == 1, TRUE,
                ifelse(ac_multiparty_g1 == 0 |
                       ac_multiparty_g2 == 0, FALSE,
                       NA)),
            ac_caliphate_direct =
                ifelse(ac_caliphate_g1 == 1 |
                       ac_caliphate_g2 == 1, TRUE,
                ifelse(ac_caliphate_g1 == 0 |
                       ac_caliphate_g2 == 0, FALSE,
                       NA)),
            ac_mandate_direct =
                ifelse(ac_mandate_g1 == 1 |
                       ac_mandate_g2 == 1, TRUE,
                ifelse(ac_mandate_g1 == 0 |
                       ac_mandate_g2 == 0, FALSE,
                       NA)),
            ## list responses to the list (groups 3 and 4; add up g1 and g2)
            ac_assad_list =
                ifelse(ac_assad_g3 >= 0, ac_assad_g3,
                ifelse(ac_assad_g4 >= 0, ac_assad_g4,
                       ac_parl_direct +
                       ac_no_parties_direct +
                       ac_mil_direct)),
            ac_caliphate_list =
                ifelse(ac_caliphate_g3 >= 0, ac_caliphate_g3,
                ifelse(ac_caliphate_g4 >= 0, ac_caliphate_g4,
                       ac_pres_direct +
                       ac_multiparty_direct +
                       ac_mandate_direct))
        ),
        ## keep only the non-helper variables
        rid,
        year,
        ac_group,
        ac_group_control,
        ac_group_treat,
        ac_group_assad,
        ac_group_caliphate,
        ac_treat_type,
        ac_assad_type,
        ac_caliphate_type,
        ac_parl_direct,
        ac_no_parties_direct,
        ac_assad_direct,
        ac_mil_direct,
        ac_pres_direct,
        ac_multiparty_direct,
        ac_caliphate_direct,
        ac_mandate_direct,
        ac_assad_list,
        ac_caliphate_list
    )

##### reassemble data

## 2017 data
AC17 <- 
    Reduce(function(...) left_join(..., by = c("rid", "year")),
           list(BG17,
                AC_list17))
## 2015 data
AC15 <- 
    Reduce(function(...) left_join(..., by = c("rid", "year")),
           list(BG15,
                AC_list15))
## combine 2015, 2017
AC <- 
    mutate(
        bind_rows(AC17, AC15),
        survey17 = ifelse(year == 2017, TRUE, FALSE)
        )

##### tinker with Lebanese province indicators

## 2017 separates out Aakar from North, so we'll recombine them to
## be consistent with 2015
## let's also pool the provinces into sensible groupings to avoid
## too-small province sample sizes
## one chunk combines the North and Aakar
## another combines Mount Lebanon and Beirut, especially since there are
## very few people in the Beirut sample
## one option combines the South and Nabatieh but leaves the Bekaa
## separate, and another combines the Bekaa with the South and Nabatieh
## on the logic that a) they're similar in practice, and b) Hizballah
## has a strong presence in all of these provinces
AC <-
    mutate(
        AC,
        ## combine North and Aakar
        province_leb_combo =
            factor(ifelse(province_leb == "Aakar",
                          which(levels(province_leb) == "North"),
                          province_leb),
                   labels =
                       levels(province_leb)[levels(province_leb) != "Aakar"]),
        ## Bekaa separate from South/Nabatieh
        province_leb_4 =
            factor(ifelse(province_leb %in% c("North", "Aakar"), 1,
                   ifelse(province_leb %in% c("Beirut", "Mount Lebanon"), 2,
                   ifelse(province_leb %in% c("South", "Nabatieh"), 3,
                   ifelse(province_leb == "Bekaa", 4, NA)))),
                   labels = c("North",
                              "Mount Lebanon",
                              "South",
                              "Bekaa")),
        ## Bekaa included in South/Nabatieh
        province_leb_3 =
            factor(ifelse(province_leb %in% c("North", "Aakar"), 1,
                   ifelse(province_leb %in% c("Beirut", "Mount Lebanon"), 2,
                   ifelse(province_leb %in% c("South", "Nabatieh", "Bekaa"),
                          3, NA))),
                   labels = c("North",
                              "Mount Lebanon",
                              "South"))
        )

#### process Arab Barometer data
##### raw data

AB4 <- read_spss("AB4-data-english.sav")

##### process data (Syrians only in Lebanon and Jordan)

ABLJ <-
    transmute(
        filter(AB4,
               as_factor(country) %in% c("Jordan", "Lebanon")),
        ## country and sample
        lebanon =
            ifelse(as_factor(country) == "Lebanon",
                   TRUE, FALSE),
        syrian  = 
            ifelse(as_factor(sample) == "Syrian Refugees",
                   TRUE, FALSE),
        province = as_factor(q1),
        district = as_factor(q2),
        ## demographics
        age     = q1001,
        older   = ifelse(age >= 30,
                         TRUE, FALSE),
        young   = !older,
        female  = ifelse(as_factor(q1002) == "Female",
                         TRUE, FALSE),
        male    = !female,
        educ    = as_factor(t1003),
        educ_hs = ifelse(educ %in% c("Secondary", "BA", "Above BA"),
                         TRUE, FALSE),
        nkids   = q1010b,
        income_low =
            ifelse(lebanon, q1015aleb, q1015ajor),
        income_low =
            ifelse(income_low > 2, NA,
            ifelse(income_low == 1,
                   TRUE, FALSE)),
        econ_status =
            factor(ifelse(q102b > 4, NA,
                          abs(q102b - 5)),
                   labels = c("Very Bad",
                              "Bad",
                              "Good",
                              "Very Good")),
        ## engagement
        understand =
            factor(ifelse(q2185 > 4, NA,
                          q2185),
                   labels = c("Low",
                              "Mid-Low",
                              "Mid-High",
                              "High")),
        interest =
            factor(ifelse(q404 > 4, NA,
                          abs(q404 - 5)),
                   labels = c("Not Interested",
                              "Somewhat Interested",
                              "Interested",
                              "Very Interested")),
        interest_bi =
            ifelse(interest %in% c("Interested",
                                   "Very Interested"),
                   TRUE, FALSE),
        ## religion 
        xian =                          # NB: only have Muslim/Christian for Syrias
            ifelse(as_factor(q1012) == "Christian",
                   TRUE, FALSE),
        religious = 
            factor(ifelse(as_factor(q609) == "Not religious", 1,
                   ifelse(as_factor(q609) == "Somewhat religious", 2,
                   ifelse(as_factor(q609) == "Religious", 3, NA))),
                   labels = c("No", "Somewhat", "Yes")),
        pray =
            factor(ifelse(q6101 > 5, NA,
                          abs(q6101 - 6)),
                   labels = c("Never",
                              "Rarely",
                              "Sometimes",
                              "Most of the time",
                              "Always")),
        quran =
            factor(ifelse(q6106 > 5, NA,
                          abs(q6106 - 6)),
                   labels = c("Never",
                              "Rarely",
                              "Sometimes",
                              "Most of the time",
                              "Always")),
        quran_high =
            ifelse(quran %in% c("Most of the time",
                                "Always"),
                   TRUE, FALSE),
        ## ISIS
        isis_threat_country =
            factor(ifelse(q826 > 3, NA,
                          abs(q826 - 4)),
                   labels = c("None",
                              "Some",
                              "Very")),
        isis_threat_arab =
            factor(ifelse(q827 > 3, NA,
                          abs(q827 - 4)),
                   labels = c("None",
                              "Some",
                              "Very")),
        isis_compatible_islam =
            factor(ifelse(q828 > 4, NA,
                          abs(q828 - 5)),
                   labels = c("Definitely No",
                              "No",
                              "Yes",
                              "Definitely Yes")),
        isis_goals = 
            factor(ifelse(q829 > 4, NA,
                          abs(q829 - 5)),
                   labels = c("Strongly Disagree",
                              "Somewhat Disagree",
                              "Somewhat Agree",
                              "Strongly Agree")),
        isis_violence = 
            factor(ifelse(q830 > 4, NA,
                          abs(q830 - 5)),
                   labels = c("Strongly Disagree",
                              "Somewhat Disagree",
                              "Somewhat Agree",
                              "Strongly Agree")),
        ## party ID
        party   =
            as_factor(q503sy),
        party_none =
            ifelse(is.na(party), NA,
            ifelse(party == "No party",
                   TRUE, FALSE)),
        party_govt =
            ifelse(is.na(party), NA,
            ifelse(party %in% c("Syrian regime",
                                "Ba'ath Party",
                                "Syrian Nationalist Party"),
                   TRUE, FALSE)),
        party_oppn =
            ifelse(party_none | party_govt, FALSE, TRUE)
        )

##### Syrians only

## filter to Syrians only, and also get rid of the empty factors
## for non-Lebanese/Jordanian provinces and districts
ABSY <-
    filter(ABLJ, syrian) %>%
    mutate(
        district = factor(district),
        province = factor(province),
        province_4 = recode(province, "El Nabatieh" = "South")
    )

### analyses
#### List experiment helper functions
##### listit code

listit <- function(Y, X, map=NULL, method="BFGS") {
  ## Returns a matrix of parameters with a row for each list item
  ## and columns for each covariate.  The parameters are placed in
  ## their appropriate cells, and the remaining cells are filled
  ## with zeros.
  place.B <- function(param, map) {
    k <- nrow(map)
    n.x <- ncol(map)
    n.param <- length(param)
    B <- matrix(NA, nrow=k,
                    ncol=n.x)
    i.param <- 1
    for (i in 1:k)
      for (j in 1:n.x) {
        if (map[i,j] == TRUE) {
          B[i,j] <- param[i.param]
          i.param <- i.param + 1
        }
        else
          B[i,j] <- 0
      }
    return(B)
  }

  neg.LL <- function(param, Y.t, Y.c, X.t, X.c, map) {
    k <- nrow(map)
    B <- place.B(param=param, map=map)
    phi <- apply((1 + exp(-X.t %*%
                          t(B[2:k,])))^(-1)
                 , 1, sum)
    p <- as.vector(((1 + exp(-X.t %*%
                             B[1,]))^(-1) +
                    phi) / k)
    q <- (1 + exp(-X.c %*% t(B[2:k,])))^(-1)
    list.LL <- sum(log(choose(k, Y.t)) +
                   Y.t * log(p) +
                   (k - Y.t) * log(1 - p))
    off.LL <- sum(Y.c * log(q) +
                  (1 - Y.c) * log(1 - q))
    joint.LL <- list.LL + off.LL
    return(-joint.LL)
  }

  neg.grad <-function(param, Y.t, Y.c, X.t, X.c, map) {
    k <- nrow(map)                      #number of list items
    n.x <- ncol(map)                    #number of covariates
    n.grad <- sum(map)                  #number of gradients to return
    B <- place.B(param=param, map=map)  #rows: list item, cols: parameters
    phi <- apply((1 + exp(-X.t %*%
                          t(B[2:k,])))^(-1)
                 , 1, sum)
    p <- as.vector(((1 + exp(-X.t %*%
                             B[1,]))^(-1) +
                    phi) / k)
    ## Calculate the derivative of p with respect to the parameters
    ## Loop through the list items
    ## For each list item, get a temporary matrix of its predictors (this.X)
    ## Use (i.grad) and (this.param) indices to make sure the results go into
    ##   the appropriate columns in the (d.p) matrix, which has a number of
    ##   rows equal to the number of treatment observations, and a number of
    ##   columns equal to the number of gradients to be calculated (i.e., the
    ##   number of parameters we're trying to estimate).
    i.grad <- 0                         #gradient index
    d.p <- matrix(0, nrow=nrow(Y.t), ncol=n.grad)
    for (i in 1:k) {
      this.X <- X.t[,1]
      this.param <- sum(map[i,])
      for (j in 2:n.x)
        if (map[i,j] == TRUE)
          this.X <- cbind(this.X, X.t[,j]) #only the Xs for this list item
      e.vec <-
        as.vector(exp(-X.t %*% B[i,]) /
                  (k * (1 + exp(-X.t %*% B[i,]))^(2)))
      d.p[,(1 + i.grad):
            (i.grad + this.param)] <-
              this.X * e.vec
      i.grad <- i.grad + this.param
    }
    ## Calculate the contributions to the gradient calculations from the
    ## off items (i.e., the non-sensitive list items).
    ## Loop through the off items (items 2 to k).
    ## For each item, get a temporary matrix of its predictors (this.X).
    ## Use (i.grad) and (this.param) indices to make sure the results go into
    ##   the appropriate column in the (off) matrix.  (Note that the (off)    
    ##   matrix leaves a number of initial columns as zeros, corresponding to
    ##   the sensitive betas which do not enter these calculations.)  The
    ##   number of rows is the number of control group observations, and the
    ##   number of columns is the number of gradients to be calculated (i.e.,
    ##   the number of parameters we're trying to estimate).
    i.grad <- sum(map[1,])              #Start after the sensitive betas
    off <- matrix(0, nrow=nrow(Y.c), ncol=n.grad)
    for (i in 2:k) {
      this.X <- X.c[,1]
      this.param <- sum(map[i,])
      for (j in 2:n.x)
        if (map[i,j] == TRUE)
          this.X <- cbind(this.X, X.c[,j])
      q <- as.vector((1 + exp(-X.c %*%
                              B[i,]))^(-1))
      off[,(1 + i.grad):(i.grad + this.param)] <-
        (Y.c[,(i - 1)] - q) * this.X
      i.grad <- i.grad + this.param
    }
    ## Calculate all the relevant components and add them together.
    ## (LL.t) is the component related to the treatment group.
    ## (LL.off) is the component related to the off items.
    ## (LL.joint) is the joint components.
    LL.t <- as.vector(Y.t * p^(-1)) * d.p +
      as.vector((k - Y.t) * (1 - p)^(-1)) * (-d.p)
    LL.t <- apply(LL.t, 2, sum)
    LL.off <- apply(off, 2, sum)
    LL.joint <- LL.t + LL.off
    return(-LL.joint)
  }

  ## Removes observations that have missing values from the
  ## working data.
  na.rm <- function(Y=Y, X=X) {
    k <- ncol(Y)                          #number of list items
    x.names <- colnames(X)                #names of the covariates
    X <- cbind(1,X)                       #append a constant to X
    colnames(X) <- c("constant", x.names) #names the columns
    ## D is a data frame.  The first column will be a vector of
    ## dummies indicating whether these observations will be used
    ## or not depending on whether or not they have NAs among the
    ## covariates.  The second column marks whether the observation
    ## is in the treatment group (ie, it gets the list) or the
    ## control group (it gets the off items individually).  The
    ## Y values follow, then the X values (with an appended constant).
    D <- as.data.frame(cbind(NA,NA,Y,X))
    colnames(D) <- c("to.use", "control",
                     colnames(Y),
                     colnames(X))
    ## If one of the off items is NA, puts NA in control, else a number
    D[,"control"] <-
      apply(D[,(4:(k + 2))], 1, sum)
    ## If control is not NA, it's a control variable, else it's NA...
    D[,"control"] <-
      ifelse(!is.na(D[,"control"]),
             TRUE, D[,"control"])
    ## ...But if there is a list count, it's not a control variable
    D[,"control"] <-
      ifelse(!is.na(D[,3]),
             FALSE, D[,"control"])
    ## Puts a count of the number of covariate NAs into to.use
    D[,"to.use"] <-
      apply(is.na(D[,((k + 3):ncol(D))]),
            1, sum)
    ## Mark to.use if there are no covariate NAs on this obs
    D[,"to.use"] <-
      ifelse((D[,"to.use"] != 0),
             FALSE, TRUE)
    ## Mark to.use as FALSE if control observations have NA Y values
    D[,"to.use"] <-
      ifelse(is.na(D[,"control"]),
             FALSE, D[,"to.use"])
    ## Take the subset of usable observations, stripping off
    ## the to.use and control columns
    D <- subset(D[3:ncol(D)], (D[,"to.use"] == TRUE))
    return(D)                             #return the usable Y and X values
  }

  ## Checks for the map defining which covariates predict which list items.
  ## If no map is supplied by the user, it assumes that every covariate in X
  ## predicts every list item.  If a map is supplied by a user, it appends
  ## an initial column of "TRUE" values to indicate that the constant is
  ## a covariate for every list item.
  if (is.null(map))
    map <- matrix(TRUE, nrow=ncol(Y),
                  ncol=(ncol(X) + 1))
  else map <- cbind(TRUE, map)
  ## Processes the data.  First it removes observations with NA
  ## values, then places into separate objects control and treatment
  ## group observations on Y and X.
  k <- ncol(Y)
  YX <- na.rm(Y=Y, X=X)
  YX.c <- subset(YX, is.na(YX[,1]))
  YX.t <- subset(YX, !is.na(YX[,1]))
  Y.c <- as.matrix(YX.c[,2:k])
  X.c <- as.matrix(YX.c[(k + 1):ncol(YX.c)])
  Y.t <- as.matrix(YX.t[,1])
  X.t <- as.matrix(YX.t[(k + 1):ncol(YX.t)])
  ## The number of parameters are the sum total TRUE values
  ## in the parameter map
  n.param <- sum(map)

  result <- optim(par=rep(0, n.param),
                  fn=neg.LL, gr=neg.grad,
                  hessian=TRUE,
                  method=method,
                  Y.t=Y.t, Y.c=Y.c,
                  X.t=X.t, X.c=X.c,
                  map=map)

  neg.se <- FALSE
  for (i in 1:ncol(X.t))
    if (is.nan(sqrt(diag(solve(result$hessian)))[i])) {
      neg.se <- TRUE
      break
    }

  output <- list(b.star=(result$par)[1:ncol(X.t)],
                 var=solve(result$hessian),
                 se=sqrt(diag(solve(result$hessian)))[1:ncol(X.t)],
                 neg.se=neg.se,
                 converged=(result$convergence == 0),
                 n.treatment=nrow(Y.t),
                 n.control=nrow(Y.c),
                 X.t=X.t,
                 items=k,
                 names=colnames(X.t),
                 logL=-result$value,
                 deviance=2 * result$value)

  class(output) <- "listit"

  return(output)
}


##### simple simulation with dummy covariates
###### simulate

## NB: x is the number of the coefficient we're working with
sim_me <- function(mod,
                   x,
                   x.seq  = 0:1,
                   p.easy = 0.10,
                   p.hard = 0.05,
                   reps   = 10) {
    ## preliminaries
    require(MASS)
    if ("listit" %in% class(mod)) {
        B <- mvrnorm(n     = reps,
                     mu    = mod$b.star,
                     Sigma = mod$var[1:length(mod$b.star),
                                     1:length(mod$b.star)])
        X <- mod$X.t
    } else {
        B <- mvrnorm(n     = reps,
                     mu    = coef(mod),
                     Sigma = vcov(mod))
        X <- model.matrix(mod)
    }
    ## simulate for each value of x
    sim <- array(data = 0, dim = c(nrow(X), reps, length(x.seq)))
    for (i in 1:length(x.seq)) {
        X[,x]    <- x.seq[i]   # set x to sequence
        sim[,,i] <- plogis(X %*% t(B))
    }
    ## return
    return(list(sim   = sim,
                x.seq = x.seq,
                reps  = reps))
}

###### calculate PPs

calc_pp <- function(sim,
                    p.easy = 0.10,
                    p.hard = 0.05) {
    ## set up
    pp <-
        matrix(data = 0,
               nrow = length(sim$x.seq),
               ncol = 5,
               dimnames = list(sim$x.seq,
                               c("lo95",
                                 "lo90",
                                 "mean",
                                 "hi90",
                                 "hi95")))
    probs <- c(p.hard/2,
               p.easy/2,
               .50,
               1 - p.easy/2,
               1 - p.hard/2)
    ## calculate predicted probabilities
    for (i in 1:length(sim$x.seq)) {
        means    <- apply(sim$sim[,,i], 2, mean)
        pp[i,]   <- quantile(means, probs = probs)
        pp[i, 3] <- mean(means)
    }
    ## return
    return(list(pp    = pp,
                sim   = sim,
                x.seq = sim$x.seq,
                reps  = sim$reps))
}

###### calculate prob differences

calc_diff <- function(sim,
                      p.easy = 0.10,
                      p.hard = 0.05) {
    ## set up
    diff <-
        matrix(data = 0,
               nrow = length(sim$x.seq) - 1,
               ncol = 5,
               dimnames = list(sim$x.seq[-1],
                               c("lo95",
                                 "lo90",
                                 "mean",
                                 "hi90",
                                 "hi95")))
    probs <- c(p.hard/2,
               p.easy/2,
               .50,
               1 - p.easy/2,
               1 - p.hard/2)
    ## calculate differences in predicted probabilities
    for (i in 1:(length(sim$x.seq) - 1)) {
        means      <- apply(sim$sim[,,i+1] -
                            sim$sim[,,i],
                            2, mean)
        diff[i,]   <- quantile(means, probs = probs)
        diff[i, 3] <- mean(means)
    }
    ## return
    return(list(diff  = diff,
                sim   = sim,
                x.seq = sim$x.seq,
                reps  = sim$reps))
}

##### simulation with two dummy covariates interacted
###### simulate interaction

sim_me_dum <- function(mod,
                       x,
                       d,
                       dx,
                       x.seq  = 0:1,
                       w.seq  = 0:1,
                       p.easy = 0.10,
                       p.hard = 0.05,
                       reps   = 10) {
    ## preliminaries
    require(MASS)
    if ("listit" %in% class(mod)) {
        B <- mvrnorm(n     = reps,
                     mu    = mod$b.star,
                     Sigma = mod$var[1:length(mod$b.star),
                                     1:length(mod$b.star)])
        X <- mod$X.t
    } else {
        B <- mvrnorm(n     = reps,
                     mu    = coef(mod),
                     Sigma = vcov(mod))
        X <- model.matrix(mod)
    }
    ## simulate for each value of x
    sim0 <- sim1 <-
        array(data = 0, dim = c(nrow(X), reps, length(x.seq)))
    for (i in 1:length(x.seq)) {
        X0 <- X                         # this data
        X0[,d] <- 1                     # turn dummies on
        X0[,c(x, dx)] <- x.seq[i]       # set x to sequence
        Xd0 <- Xd1 <- X0                # dummied datasets
        Xd0[,c(d, dx)] <- 0             # turn off dummies
        sim0[,,i] <- plogis(Xd0 %*% t(B))
        sim1[,,i] <- plogis(Xd1 %*% t(B))
    }
    ## return
    return(list(sim0   = sim0,
                sim1   = sim1,
                x.seq = x.seq,
                reps  = reps))
}

###### calculate PPs interaction

calc_pp_dum <- function(sim,
                        p.easy = 0.10,
                        p.hard = 0.05) {
    ## set up
    pp0 <- pp1 <-
        matrix(data = 0,
               nrow = length(sim$x.seq),
               ncol = 5,
               dimnames = list(sim$x.seq,
                               c("lo95",
                                 "lo90",
                                 "mean",
                                 "hi90",
                                 "hi95")))
    probs <- c(p.hard/2,
               p.easy/2,
               .50,
               1 - p.easy/2,
               1 - p.hard/2)
    ## calculate predicted probabilities
    for (i in 1:length(sim$x.seq)) {
        means0    <- apply(sim$sim0[,,i], 2, mean)
        pp0[i,]   <- quantile(means0, probs = probs)
        pp0[i, 3] <- mean(means0)
        means1    <- apply(sim$sim1[,,i], 2, mean)
        pp1[i,]   <- quantile(means1, probs = probs)
        pp1[i, 3] <- mean(means1)
    }
    ## return
    return(list(pp0   = pp0,
                pp1   = pp1,
                sim   = sim,
                x.seq = sim$x.seq,
                reps  = sim$reps))
}

###### calculate prob differences interaction

calc_diff_dum <- function(sim,
                          p.easy = 0.10,
                          p.hard = 0.05) {
    ## set up
    diff0 <- diff1 <-
        matrix(data = 0,
               nrow = length(sim$x.seq) - 1,
               ncol = 5,
               dimnames = list(sim$x.seq[-1],
                               c("lo95",
                                 "lo90",
                                 "mean",
                                 "hi90",
                                 "hi95")))
    probs <- c(p.hard/2,
               p.easy/2,
               .50,
               1 - p.easy/2,
               1 - p.hard/2)
    ## calculate differences in predicted probabilities
    for (i in 1:(length(sim$x.seq) - 1)) {
        means0      <- apply(sim$sim0[,,i+1] -
                             sim$sim0[,,i],
                             2, mean)
        diff0[i,]   <- quantile(means0, probs = probs)
        diff0[i, 3] <- mean(means0)
        means1      <- apply(sim$sim1[,,i+1] -
                             sim$sim1[,,i],
                             2, mean)
        diff1[i,]   <- quantile(means1, probs = probs)
        diff1[i, 3] <- mean(means1)
    }
    ## return
    return(list(diff0 = diff0,
                diff1 = diff1,
                sim   = sim,
                x.seq = sim$x.seq,
                reps  = sim$reps))
}

#### List experiment analyses
##### Assad
###### set up container list

## gathers all the result tables in a tidy list
## base: distinguish good/okay, 2015/2017
## pool: pool good/okay, 2015/2017
## pall: pool all: pool good/okay, pool 2015/2017
## example access:
##  - Assad$pool$main$full
##    - pooled good/okay metric
##    - main direct/list results
##    - full sample
##  - Assad$pool$diff$full
##    - pooled good/okay metric
##    - difference between the list and direct results
##    - full sample
Assad <-
    list(base = list(main = list(),
                     diff = list()),
         pool = list(main = list(),
                     diff = list()),
         pall = list(main = list(),
                     diff = list()))

###### good/okay, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = ac_assad_type,
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$base$main$full <- Both
Assad$base$diff$full <- Diff

###### good/okay, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = ac_assad_type,
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$base$main$sub0 <- Both
Assad$base$diff$sub0 <- Diff

###### good/okay, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = ac_assad_type,
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$base$main$sub1 <- Both
Assad$base$diff$sub1 <- Diff

###### pooled, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pool$main$full <- Both
Assad$pool$diff$full <- Diff

###### pooled, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pool$main$sub0 <- Both
Assad$pool$diff$sub0 <- Diff

###### pooled, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pool$main$sub1 <- Both
Assad$pool$diff$sub1 <- Diff

###### pool metric and year, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pall$main$full <- Both
Assad$pall$diff$full <- Diff

###### pool metric and year, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pall$main$sub0 <- Both
Assad$pall$diff$sub0 <- Diff

###### pool metric and year, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_assad_direct,
                    q_list    = ac_assad_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Assad$pall$main$sub1 <- Both
Assad$pall$diff$sub1 <- Diff

##### Jihad
###### set up container list

## gathers all the result tables in a tidy list
## base: distinguish good/okay, 2015/2017
## pool: pool good/okay, 2015/2017
## pall: pool all: pool good/okay, pool 2015/2017
## example access:
##  - Jihad$pool$main$full
##    - pooled good/okay metric
##    - main direct/list results
##    - full sample
##  - Jihad$pool$diff$full
##    - pooled good/okay metric
##    - difference between the list and direct results
##    - full sample
Jihad <-
    list(base = list(main = list(),
                     diff = list()),
         pool = list(main = list(),
                     diff = list()),
         pall = list(main = list(),
                     diff = list()))

###### good/okay, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = ac_caliphate_type,
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$base$main$full <- Both
Jihad$base$diff$full <- Diff

###### good/okay, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = ac_caliphate_type,
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$base$main$sub0 <- Both
Jihad$base$diff$sub0 <- Diff

###### good/okay, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = ac_caliphate_type,
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$base$main$sub1 <- Both
Jihad$base$diff$sub1 <- Diff

###### pooled, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pool$main$full <- Both
Jihad$pool$diff$full <- Diff

###### pooled, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pool$main$sub0 <- Both
Jihad$pool$diff$sub0 <- Diff

###### pooled, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pool$main$sub1 <- Both
Jihad$pool$diff$sub1 <- Diff

###### pool metric and year, full sample
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(1, labels = "full")),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pall$main$full <- Both
Jihad$pall$diff$full <- Diff

###### pool metric and year, govt/oppn/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_4cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pall$main$sub0 <- Both
Jihad$pall$diff$sub0 <- Diff

###### pool metric and year, govt/fsa/islamists/none
####### calculate

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = fac_5cat),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pall$main$sub1 <- Both
Jihad$pall$diff$sub1 <- Diff

###### pool metric and year, govt and none/oppn
####### calculate

table(AC$fac_4cat)

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(ifelse(fac_4cat == "oppn",
                                              1, 2),
                                       labels = c("oppn",
                                                  "non-oppn"))),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

####### results

Jihad$pall$main$sub2 <- Both
Jihad$pall$diff$sub2 <- Diff

##### Ad Hoc
###### Assad among government sympathizers by private religion

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = factor(quran_weekly, labels = c("Low", "High"))
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

###### Assad among government sympathizers by public religion

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = factor(relpol_bi, labels = c("Low", "High"))
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

###### Assad among government sympathizers by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

###### Assad among opposition sympathizers by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_oppn) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )
###### Assad among unaligned by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_none) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )
###### Jihadis among FSA/domestic/foreign pool metric and year

## gather data
Dat0 <-
    group_by(mutate(AC,
                    q_direct  = ac_caliphate_direct,
                    q_list    = ac_caliphate_list,
                    metric    = factor(1, labels = "pool"),
                    year      = factor(1, labels = "pool"),
                    treat     = ac_treat_type,
                    faction   = factor(ifelse(fac_all == "Nationalists", 1,
                                       ifelse(fac_all == "Domestic Islamists", 2,
                                       ifelse(fac_all == "Foreign Islamists", 3,
                                              NA))),
                                       labels = c("Nationalists",
                                                  "Domestic Islamists",
                                                  "Foreign Islamists"))),
             year,
             metric,
             treat,
             faction)
## Dat0

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Reblist <- Both[5:7, c("faction", "mean", "se")]

Reblist

## foreign islamists vs nationalists
m  <- unlist(Reblist[3, 2] - Reblist[1, 2])           # diff in means
se <- unlist(sqrt(Reblist[3, 3]^2 + Reblist[1, 3]^2)) # SE(DIM)
c(mean = m,
  se   = se,
  p    = pnorm(-m/se) * 2)

## domestic islamists vs nationalists
m  <- unlist(Reblist[2, 2] - Reblist[1, 2])           # diff in means
se <- unlist(sqrt(Reblist[2, 3]^2 + Reblist[1, 3]^2)) # SE(DIM)
c(mean = m,
  se   = se,
  p    = pnorm(-m/se) * 2)

## foreign islamists vs domestic islamists
m  <- unlist(Reblist[3, 2] - Reblist[2, 2])           # diff in means
se <- unlist(sqrt(Reblist[3, 3]^2 + Reblist[2, 3]^2)) # SE(DIM)
c(mean = m,
  se   = se,
  p    = pnorm(-m/se) * 2)

###### Jihadis among government sympathizers by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_caliphate_direct,
        q_list    = ac_caliphate_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

###### Jihadis among opposition sympathizers by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_oppn) %>%
    mutate(
        q_direct  = ac_caliphate_direct,
        q_list    = ac_caliphate_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

###### Jihadis among unaligned by Lebanese province

## gather data
Dat0 <-
    AC %>%
    filter(fac_none) %>%
    mutate(
        q_direct  = ac_caliphate_direct,
        q_list    = ac_caliphate_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = province_leb_3
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Ctrl

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )
## Treat

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]
## Direct

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)
## Indirect

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))
## Both

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))
## Diff

Both %>%
    filter(!is.na(faction))

Diff %>%
    filter(!is.na(faction))

Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )
### main text
#### Table 1: summary statistics table
##### Syria pooled

## gather and tidy summary statistics for the Syria data,
## pooling 2015 and 2017 together

## I know dplyr is great and all, but it sure is inscrutable at times:
## https://stackoverflow.com/questions/34594641/dplyr-summary-table-for-multiple-variables

sumstats <-
    summarize_all(
        select(
            mutate(AC,
                   income = NA,
                   young  = age < 30,
                   srd2p  = syr_room_density >= 2,
                   lrd2p  = leb_room_density >= 2),
            income,
            young,
            male,
            minority,
            srd2p,
            lrd2p,
            educ     = educ_hs,
            interest = interest_bi,
            quran    = quran_weekly,
            relpol   = relpol_bi
        ),
        funs(mean = mean(., na.rm = TRUE),
             sd   = sd(., na.rm = TRUE),
             n    = sum(!is.na(.)))
    )

sstidy_pool  <-
    sumstats %>%
    gather(stat, val) %>%
    separate(stat, into = c("var", "stat"), sep = "_") %>%
    spread(stat, val) %>%
    select(var, mean, sd, n)

##### Syria 2015

## gather and tidy summary statistics for the 2015 Syria data

## I know dplyr is great and all, but it sure is inscrutable at times:
## https://stackoverflow.com/questions/34594641/dplyr-summary-table-for-multiple-variables

sumstats <-
    summarize_all(
        select(
            mutate(filter(AC, year == 2015),
                   income = NA,
                   young  = age < 30,
                   srd2p  = syr_room_density >= 2,
                   lrd2p  = leb_room_density >= 2),
            income,
            young,
            male,
            minority,
            srd2p,
            lrd2p,
            educ     = educ_hs,
            interest = interest_bi,
            quran    = quran_weekly,
            relpol   = relpol_bi
        ),
        funs(mean = mean(., na.rm = TRUE),
             sd   = sd(., na.rm = TRUE),
             n    = sum(!is.na(.)))
    )

sstidy_2015  <-
    sumstats %>%
    gather(stat, val) %>%
    separate(stat, into = c("var", "stat"), sep = "_") %>%
    spread(stat, val) %>%
    select(var, mean, sd, n)

##### Syria 2017

## gather and tidy summary statistics for the 2017 Syria data

## I know dplyr is great and all, but it sure is inscrutable at times:
## https://stackoverflow.com/questions/34594641/dplyr-summary-table-for-multiple-variables

sumstats <-
    summarize_all(
        select(
            mutate(filter(AC, year == 2017),
                   income = NA,
                   young  = age < 30,
                   srd2p  = syr_room_density >= 2,
                   lrd2p  = leb_room_density >= 2),
            income,
            young,
            male,
            minority,
            srd2p,
            lrd2p,
            educ     = educ_hs,
            interest = interest_bi,
            quran    = quran_weekly,
            relpol   = relpol_bi
        ),
        funs(mean = mean(., na.rm = TRUE),
             sd   = sd(., na.rm = TRUE),
             n    = sum(!is.na(.)))
    )

sstidy_2017  <-
    sumstats %>%
    gather(stat, val) %>%
    separate(stat, into = c("var", "stat"), sep = "_") %>%
    spread(stat, val) %>%
    select(var, mean, sd, n)

##### AB Lebanon

## gather and tidy summary statistics for the Lebanon Syrian data
## from the Arab Barometer

## I know dplyr is great and all, but it sure is inscrutable at times:
## https://stackoverflow.com/questions/34594641/dplyr-summary-table-for-multiple-variables

sumstats <-
    summarize_all(
        select(
            mutate(filter(ABSY, lebanon),
                   srd2p    = NA,
                   lrd2p    = NA,
                   relpol   = NA,
                   ),
            srd2p,
            lrd2p,
            relpol,
            young,
            male,
            minority = xian,
            income   = income_low,
            educ     = educ_hs,
            interest = interest_bi,
            quran    = quran_high
        ),
        funs(mean = mean(., na.rm = TRUE),
             sd   = sd(., na.rm = TRUE),
             n    = sum(!is.na(.)))
    )

sstidy_leb <-
    sumstats %>%
    gather(stat, val) %>%
    separate(stat, into = c("var", "stat"), sep = "_") %>%
    spread(stat, val) %>%
    select(var, mean, sd, n)

##### AB Jordan

## gather and tidy summary statistics for the Jordan Syrian data
## from the Arab Barometer

## I know dplyr is great and all, but it sure is inscrutable at times:
## https://stackoverflow.com/questions/34594641/dplyr-summary-table-for-multiple-variables

sumstats <-
    summarize_all(
        select(
            mutate(filter(ABSY, !lebanon),
                   srd2p    = NA,
                   lrd2p    = NA,
                   relpol   = NA,
                   ),
            srd2p,
            lrd2p,
            relpol,
            young,
            male,
            minority = xian,
            income   = income_low,
            educ     = educ_hs,
            interest = interest_bi,
            quran    = quran_high
        ),
        funs(mean = mean(., na.rm = TRUE),
             sd   = sd(., na.rm = TRUE),
             n    = sum(!is.na(.)))
    )

sstidy_jor <-
    sumstats %>%
    gather(stat, val) %>%
    separate(stat, into = c("var", "stat"), sep = "_") %>%
    spread(stat, val) %>%
    select(var, mean, sd, n)

##### gather

## now we need to smoosh everything together into a data frame,
## renooberate it with dplyr, and make some pretty labels

sstidy <-
    mutate(bind_rows(sstidy_pool,
                     sstidy_2015,
                     sstidy_2017,
                     sstidy_leb,
                     sstidy_jor),
           sample =
               rep(factor(1:5,
                          labels = c("pool",
                                     "2015",
                                     "2017",
                                     "leb",
                                     "jor")),
                   each = n() / 5))

sumstats <- 
    select(sstidy,
           var,
           sample,
           mean,
           sd,
           n) %>%
    gather(varname,
           varval,
           mean:n) %>%
    unite(tempvar,
          varname,
          sample) %>%
    spread(tempvar,
           varval) %>%
    select(var,
           mean_pool,
           sd_pool,
           n_pool,
           mean_2015,
           sd_2015,
           n_2015,
           mean_2017,
           sd_2017,
           n_2017,
           mean_leb,
           sd_leb,
           n_leb,
           mean_jor,
           sd_jor,
           n_jor) %>%
    arrange(match(var,
                  c("young",
                    "male",
                    "minority",
                    "srd2p",
                    "lrd2p",
                    "income",
                    "educ",
                    "interest",
                    "quran",
                    "relpol")))

## descriptive variable labels
sumstats$var <-
    c("Young (Age $< 30$)",
      "Male",
      "Minority (Kurds, Non-Sunnis)",
      "Crowded in Syria (Room Density $\\ge 2$)",
      "Crowded in Lebanon (Room Density $\\ge 2$)",
      "Low Income",
      "Educated (Completed Secondary School)",
      "Interested in Politics",
      "Private Religion (Reads Quran/Bible Weekly)",
      "Public Religion (Role for Religion in Public Affairs)")

## only want the means and mean differences
## (they're dummies, so the SD is implied,
##  and NAs are trivial in number)
sumstats_useme <-
    select(sumstats,
           var,
           mean_pool,
           mean_2015,
           mean_2017,
           mean_leb,
           mean_jor)

##### print table to LaTeX

## and here we go: let's print out everything into a pretty LaTeX table
## it'll print out to the R console, and we can drop it right into
## the .tex file

## convert proportions to percentages
sumstats_useme[,-1] <- sumstats_useme[,-1] * 100

## print it out
print(xtable(sumstats_useme,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 5)),
             digits  = 1,
             caption = paste("Background Summary Statistics ",
                             "(Percentages from Dichotomous Variables)",
                             sep = "\n"),
             label   = "tab:sumstats"),
      table.placement  = "htbp",
      ## ugh: getting a pretty table header is a mess
      add.to.row       = list(pos     = list(0),
                              command =
                                  ## samples, midrules, sumstat names
                                  paste0(
                                      paste0("& ",
                                             paste("\\multicolumn{3}{c}{Original Data}",
                                                   "\\multicolumn{2}{c}{Arab Barometer}",
                                                   sep = " & "),
                                             " \\\\ \n",
                                             collapse = ""),
                                      paste0("\\cmidrule(lr){2-4}",
                                             "\\cmidrule(lr){5-6}",
                                             "\n",
                                             collapse = " "),
                                      paste0("& ", # samples
                                             paste("\\multicolumn{1}{c}{Pooled}",
                                                   "\\multicolumn{1}{c}{2015}",
                                                   "\\multicolumn{1}{c}{2017}",
                                                   "\\multicolumn{1}{c}{Lebanon}",
                                                   "\\multicolumn{1}{c}{Jordan}",
                                                   sep = " & "),
                                             " \\\\ \n",
                                             collapse = "")
                                  )),
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

#### Table 2: declared support for factions

## fiddle with data/labels, drop "other"
Dat0 <-
    mutate(AC,
           fac3 =
               factor(ifelse(fac_4cat == "none", 1,
                      ifelse(fac_4cat == "govt", 2,
                      ifelse(fac_4cat == "oppn", 3,
                             NA))),
                      labels = c("Unaligned",
                                 "Government",
                                 "Rebels")),
           fac4 =
               factor(ifelse(fac_5cat == "none", 1,
                      ifelse(fac_5cat == "govt", 2,
                      ifelse(fac_5cat == "fsa", 3,
                      ifelse(fac_5cat == "islamists", 4,
                             NA)))),
                      labels = c("Unaligned",
                                 "Government",
                                 "Nationalists",
                                 "Islamists")),
           fac5 =
               factor(ifelse(fac_all == "No Preference", 1,
                      ifelse(fac_all == "Government", 2,
                      ifelse(fac_all == "Nationalists", 3,
                      ifelse(fac_all == "Domestic Islamists", 4,
                      ifelse(fac_all == "Foreign Islamists", 5,
                             NA))))),
                      labels = c("Unaligned",
                                 "Government",
                                 "Nationalists",
                                 "Domestic Islamists",
                                 "Foreign Islamists")))

## unaligned, government, and rebels
factab3 <- 
    data_frame(Faction =
                   levels(Dat0$fac3),
               Pooled =
                   as.vector(prop.table(table(Dat0$fac3))),
               year2015 =
                   with(filter(Dat0, year == 2015),
                        as.vector(prop.table(table(fac3)))),
               year2017 =
                   with(filter(Dat0, year == 2017),
                        as.vector(prop.table(table(fac3)))),
               diff =
                   year2017 - year2015)
## separate out different kinds of rebels
factab4 <- 
    data_frame(Faction =
                   levels(Dat0$fac4),
               Pooled =
                   as.vector(prop.table(table(Dat0$fac4))),
               year2015 =
                   with(filter(Dat0, year == 2015),
                        as.vector(prop.table(table(fac4)))),
               year2017 =
                   with(filter(Dat0, year == 2017),
                        as.vector(prop.table(table(fac4)))),
               diff =
                   year2017 - year2015)
## separate out different kinds of Islamists
factab5 <- 
    data_frame(Faction =
                   levels(Dat0$fac5),
               Pooled =
                   as.vector(prop.table(table(Dat0$fac5))),
               year2015 =
                   with(filter(Dat0, year == 2015),
                        as.vector(prop.table(table(fac5)))),
               year2017 =
                   with(filter(Dat0, year == 2017),
                        as.vector(prop.table(table(fac5)))),
               diff =
                   year2017 - year2015)

## bind the relevant data
factab0 <-
    bind_rows(factab3,
              factab4[3:4,],
              factab5[4:5,])
factab0$Faction <-
    with(factab0,
         ifelse(Faction %in% c("Nationalists",
                               "Islamists"),
                paste("$\\cdots$", Faction),
         ifelse(Faction %in% c("Domestic Islamists",
                               "Foreign Islamists"),
                paste("$\\cdots\\cdots$", Faction),
                Faction)))
## convert to percentages
factab0[,-1] <- factab0[,-1] * 100

## print it out (simple-ish)
print(setNames(xtable(factab0,
                      align   = c("l ",
                                  "l ", # extra for (removed) row numbers
                                  rep("D{.}{.}{2.2} ", 4)),
                      digits  = 1,
                      caption = "Percent Support for Belligerents Derived from Rank Ordering",
                      label   = "tab:rank"),
               c("Faction",
                 "Pooled",
                 "2015",
                 "2017",
                 "Change")),
      type                   = "latex",
      table.placement        = "htbp",
      sanitize.text.function = function(x){x},
      ## ugh: getting a pretty table header is a mess
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # column headers
                                         paste0("\\multicolumn{1}{c}{",
                                                c("Pooled", "2015", "2017", "Change"),
                                                "}",
                                                collapse = " & "),
                                         " \\\\ \n")),
      include.rownames = FALSE,
      include.colnames = FALSE,
      booktabs         = TRUE)

## print it out (fancy!)
print(setNames(xtable(factab0,
                      align   = c("l ",
                                  "l ", # extra for (removed) row numbers
                                  rep("D{.}{.}{2.2} ", 4)),
                      digits  = 1,
                      caption = "Percent Declared Support for Belligerents Derived from Rank Ordering",
                      label   = "tab:rank"),
               c("Faction",
                 "Pooled",
                 "2015",
                 "2017",
                 "Change")),
      type                   = "latex",
      table.placement        = "htbp",
      sanitize.text.function = function(x){x},
      ## ugh: getting a pretty table header is a mess
      ## prettify with a horizontal line and extra "Rebel" header
      ## to break up the Unaligned/Govt/Oppn chunk from the oppn subgroups
      add.to.row       = list(pos     = list(0, 3),
                              command =
                                  c(paste0("& ", # column headers
                                           paste0("\\multicolumn{1}{c}{",
                                                  c("Pooled", "2015", "2017", "Change"),
                                                  "}",
                                                  collapse = " & "),
                                           " \\\\ \n"),
                                    paste0("\\midrule\n",
                                           "\\multicolumn{5}{l}{Rebels} \\\\ \n"))),
      include.rownames = FALSE,
      include.colnames = FALSE,
      booktabs         = TRUE)

#### Figure 1: Assad and Jihadis, pooled samples and metrics

## gather list experiment data
Dat0 <-
    bind_rows(
        ## Assad data
        transmute(
            Assad$pall$main$full,
            mean,
            se,
            dv = factor(1, labels = "Assad"),
            ## put direct Q above indirect Q in plot
            qtype = factor(qtype,
                           levels = levels(qtype)[2:1],
                           labels = c("List",
                                      "Direct"))
        ),
        ## Jihadi data
        transmute(
            Jihad$pall$main$full,
            mean,
            se,
            dv = factor(2, labels = "Jihadis"),
            ## put direct Q above indirect Q in plot
            qtype = factor(qtype,
                           levels = levels(qtype)[2:1],
                           labels = c("List",
                                      "Direct"))
        )
    )

## plot it
plot_both_poolall <-
    ggplot(Dat0,
           aes(x = qtype,
               y = mean)) +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-.05, .55),
                       labels       = seq(0.0, 1.0, .25) * 100,
                       breaks       = seq(0.0, 1.0, .25),
                       minor_breaks = seq(-0.125, 1.125, .25)) +
    ## misc
    coord_flip() +
    facet_wrap(~ dv, ncol = 2) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_both_poolall_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_both_poolall),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_both_poolall_egg,
       filename = "both-poolall.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 1.50)

#### Figure 2: Assad by faction, pooled samples and metrics

## gather list experiment data
Dat0 <-
    transmute(
        filter(Assad$pall$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        ## faction = as.character(faction),
        ## dv = factor(1, labels = "Assad"),
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct")),
        ## highlight facets that have direct/indirect differences
        high    = ifelse(faction == "Rebels",
                         FALSE,
                         TRUE)
    )

## plot it
plot_assad_poolall <-
    ggplot(Dat0,
           aes(x    = qtype,
               y    = mean,
               fill = high)) +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## highlight facets
    geom_rect(xmin  = -Inf,
              xmax  =  Inf,
              ymin  = -Inf,
              ymax  =  Inf,
              alpha = 0.035) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.10, 1.10),
                       labels       = seq(0.0, 1.0, .50) * 100,
                       breaks       = seq(0.0, 1.0, .50),
                       minor_breaks = seq(0.25, 0.75, .25)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_wrap(~ faction, nrow = 1) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_assad_poolall_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_assad_poolall),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_assad_poolall_egg,
       filename = "assad-poolall.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 1.50)

#### Figure 3: Jihadis by faction, pooled samples and metrics

## gather list experiment data
Dat0 <-
    transmute(
        filter(Jihad$pall$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        ## faction = as.character(faction),
        ## dv = factor(1, labels = "Assad"),
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct")),
        ## highlight estimates with negative probabilities
        neg     = ifelse(mean < 0 &
                         pnorm(-abs(mean/se)) * 2 < .15,
                         TRUE, FALSE),
        ## highlight facets that have direct/indirect differences
        high    = ifelse(faction == "Rebels",
                         FALSE,
                         TRUE)
    )

## plot it
plot_jihad_poolall <-
    ggplot(Dat0,
           aes(x     = qtype,
               y     = mean,
               color = neg,
               fill  = high)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## highlight facets
    geom_rect(xmin  = -Inf,
              xmax  =  Inf,
              ymin  = -Inf,
              ymax  =  Inf,
              alpha = 0.035) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.27, 0.27),
                       labels       = seq(-1.0, 1.0, .20) * 100,
                       breaks       = seq(-1.0, 1.0, .20),
                       minor_breaks = seq(-1.0, 1.0, .10)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    scale_color_manual(values = c("black", "red4")) +
    coord_flip() +
    facet_wrap(~ faction, nrow = 1) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_jihad_poolall_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_jihad_poolall),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_jihad_poolall_egg,
       filename = "jihadi-poolall-3cat.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 1.50)

#### Figure 4: Assad change in support over time by faction
##### gather list data

## gather list experiment data
Dat_list <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               qtype   == "list",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )

Dat_direct <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               qtype   == "direct",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )


## ugh, break this up
L_Govt <- subset(Dat_list, faction == "Government")
L_Rebs <- subset(Dat_list, faction == "Rebels")
L_None <- subset(Dat_list, faction == "Unaligned")

D_Govt <- subset(Dat_direct, faction == "Government")
D_Rebs <- subset(Dat_direct, faction == "Rebels")
D_None <- subset(Dat_direct, faction == "Unaligned")

L_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(L_Govt$mean * c(-1, 1)),
                     sum(L_Rebs$mean * c(-1, 1)),
                     sum(L_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(L_Govt$se^2)),
                     sqrt(sum(L_Rebs$se^2)),
                     sqrt(sum(L_None$se^2))))

D_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(D_Govt$mean * c(-1, 1)),
                     sum(D_Rebs$mean * c(-1, 1)),
                     sum(D_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(D_Govt$se^2)),
                     sqrt(sum(D_Rebs$se^2)),
                     sqrt(sum(D_None$se^2))))

DJ_Diff <-
    mutate(
        bind_rows(D_Diff,
                  L_Diff),
        dv =
            factor(rep(2:1, each = 3),
                   labels = c("List",
                              "Direct")),
        high_point =
            ifelse((faction == "Government" & dv == "List") |
                   (faction == "Unaligned"  & dv == "Direct"),
                   TRUE,
                   FALSE),
        high =
            ifelse(faction == "Rebels",
                   FALSE,
                   TRUE)
    )

##### plot it

## plot it
plot_assad_yeardiff <-
    ggplot(DJ_Diff,
           aes(x     = dv,
               y     = mean,
               color = high_point,
               fill  = high)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## highlight facets
    geom_rect(xmin  = -Inf,
              xmax  =  Inf,
              ymin  = -Inf,
              ymax  =  Inf,
              alpha = 0.035) +
    ## percent scales
    labs(x = NULL,
         y = "Change in Percent Support") +
    scale_y_continuous(limits       = c(-0.40, 0.40),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-0.875, 0.875, .125)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    scale_color_manual(values = c("black", "red4")) +
    coord_flip() +
    facet_wrap(~ faction, ncol = 3) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_assad_yeardiff_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_assad_yeardiff),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_assad_yeardiff_egg,
       filename = "assad-yeardiff-3cat.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 1.50)

#### Figure 5: Increase in pro-Assad dissimulation
##### gather data

Dat_Assad <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct"))
    )

##### process data

## we want the *change* in misrepresentation over surveys:
## (List15 - Direct15) - (List17 - Direct17)
## that's where the c(-1, 1, 1, -1) business comes from below

## ugh, break this up
A_Govt <- subset(Dat_Assad, faction == "Government")
A_Rebs <- subset(Dat_Assad, faction == "Rebels")
A_None <- subset(Dat_Assad, faction == "Unaligned")

A_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(A_Govt$mean * c(-1, 1, 1, -1)),
                     sum(A_Rebs$mean * c(-1, 1, 1, -1)),
                     sum(A_None$mean * c(-1, 1, 1, -1))),
               se =
                   c(sqrt(sum(A_Govt$se^2)),
                     sqrt(sum(A_Rebs$se^2)),
                     sqrt(sum(A_None$se^2))),
               high =
                   ifelse(faction == "Rebels",
                          FALSE,
                          TRUE))

## reorder factors according to means from Assad outcome
A_Diff$faction <-
    reorder(A_Diff$faction, -A_Diff$mean)

##### plot it

## plot it
plot_misrep_assad <-
    ggplot(A_Diff,
           aes(x     = faction,
               y     = mean,
               color = high)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Change in Direct - List") +
    scale_y_continuous(limits       = c(-0.15, 0.90),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-0.875, 0.875, .125)) +
    ## misc
    scale_color_manual(values = c("black", "red4")) +
    coord_flip() +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_misrep_assad_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_misrep_assad),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_misrep_assad_egg,
       filename = "misrep-assad.pdf",
       family   = "Palatino",
       width    = 2.50,
       height   = 1.00)

#### Footnote 1: design effects footnote
##### preliminaries

## ict.test is fussy, so we need to omit NAs

## Assad
ICT_assad_1517 <-
    na.omit(transmute(AC,
                      assad = ac_assad_list,
                      treat = as.numeric(ac_group_treat)))
ICT_assad_15 <-
    na.omit(transmute(filter(AC, year == 2015),
                      assad = ac_assad_list,
                      treat = as.numeric(ac_group_treat)))
ICT_assad_17 <-
    na.omit(transmute(filter(AC, year == 2017),
                      assad = ac_assad_list,
                      treat = as.numeric(ac_group_treat)))

## jihadis
ICT_jihad_1517 <-
    na.omit(transmute(AC,
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))
ICT_jihad_15 <-
    na.omit(transmute(filter(AC, year == 2015),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))
ICT_jihad_17 <-
    na.omit(transmute(filter(AC, year == 2017),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))

## let's break down the jihadis by faction and year
ICT_jihad_15_oppn <-
    na.omit(transmute(filter(AC,
                             fac_4cat == "oppn",
                             year == 2015),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))
ICT_jihad_15_nonoppn <-
    na.omit(transmute(filter(AC,
                             fac_4cat != "oppn",
                             year == 2015),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))
ICT_jihad_17_oppn <-
    na.omit(transmute(filter(AC,
                             fac_4cat == "oppn",
                             year == 2017),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))
ICT_jihad_17_nonoppn <-
    na.omit(transmute(filter(AC,
                             fac_4cat != "oppn",
                             year == 2017),
                      jihad = ac_caliphate_list,
                      treat = as.numeric(ac_group_treat)))

##### design effect tests

## do the test on the Assad and jihadi lists in the pooled sample,
## in 2015, and in 2017, and display the p-values of the test
## shows no design effects for Assad (p >= .987), but shows a
## design effect for the jihadis in 2017 (p = .0191)
Design <-
    data_frame(
        outcome =
            factor(rep(1:2, each = 3),
                   labels = c("Assad", "Jihad")),
        sample =
            factor(rep(1:3, 2),
                   labels = c("Both", "2015", "2017")),
        pval =
            c(with(ICT_assad_1517, ict.test(assad, treat, 3))$p,
              with(ICT_assad_15,   ict.test(assad, treat, 3))$p,
              with(ICT_assad_17,   ict.test(assad, treat, 3))$p,
              with(ICT_jihad_1517, ict.test(jihad, treat, 3))$p,
              with(ICT_jihad_15,   ict.test(jihad, treat, 3))$p,
              with(ICT_jihad_17,   ict.test(jihad, treat, 3))$p)
    )
Design

## do the test for the jihadi list in 2015 and 2017 among rebel and
## non-rebel subsamples
## shows a design effect among non-rebels in both years (p < .005)
Design_jihad <- 
    data_frame(
        faction =
            factor(rep(1:2, 2),
                   labels = c("rebels", "non-rebels")),
        sample =
            factor(rep(1:2, each = 2),
                   labels = c("2015", "2017")),
        pval =
            c(with(ICT_jihad_15_oppn,    ict.test(jihad, treat, 3))$p,
              with(ICT_jihad_15_nonoppn, ict.test(jihad, treat, 3))$p, 
              with(ICT_jihad_17_oppn,    ict.test(jihad, treat, 3))$p, 
              with(ICT_jihad_17_nonoppn, ict.test(jihad, treat, 3))$p)
        )
Design_jihad

#### Footnote 4: religious support for Assad footnote
##### Assad among government sympathizers by private religion

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = factor(quran_weekly, labels = c("Low", "High"))
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))

## print out the results:
## "faction" here is actually private religion (low vs high)
## print out the mean and the 95% plus/minus values, along with
## the 95% confidence bounds
## top two rows are for the direct question, bottom two rows
## are for the indirect list question
## top two rows show nearly 100 percent approval regardless
## of level of religiosity
## bottom two rows show 48.5 percent approval for low religiosity,
## and 28.5 percent approval for high religiosity
Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

##### Assad among government sympathizers by public religion

## gather data
Dat0 <-
    AC %>%
    filter(fac_govt) %>%
    mutate(
        q_direct  = ac_assad_direct,
        q_list    = ac_assad_list,
        metric    = factor(1, labels = "pool"),
        year      = factor(1, labels = "pool"),
        treat     = ac_treat_type,
        faction   = factor(relpol_bi, labels = c("Low", "High"))
    ) %>%
    group_by(
        year,
        metric,
        treat,
        faction
    )

## raw control group summary
Ctrl <-
    summarize(
        filter(Dat0, treat == "ctrl"),
        mean  = mean(q_direct, na.rm = TRUE),
        sd    = sqrt(mean * (1 - mean)),
        n     = n(),
        n_use = sum(!is.na(q_direct)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )

## raw treatment group summary
Treat <- 
    summarize(
        Dat0,
        mean  = mean(q_list, na.rm = TRUE),
        sd    = sd(q_list, na.rm = TRUE),
        n     = n(),
        n_use = sum(!is.na(q_list)),
        na    = n - n_use,
        na_pr = na / n,
        se    = sd / sqrt(n_use)
    )

## direct Q summary
## NOTE (grumble grumble):
##  - dplyr is going to add "treat" if I don't
##  - remove the "treat" column with the [,-1]
Direct  <-
    select(
        Ctrl,
        treat,
        year,
        metric,
        faction,
        mean,
        se,
        n,
        n_use,
        na_pr
    )[,-1]

## list Q estimate and SE
Indirect <- 
    select(Treat,
           year,
           metric,
           treat,
           faction,
           mean,
           se,
           n,
           n_use,
           na_pr) %>%
    gather(varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          treat,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(faction,
              mean   = treat_mean - ctrl_mean,
              se     = sqrt(treat_se^2 + ctrl_se^2),
              n      = treat_n,
              n_use  = treat_n_use,
              na_pr  = treat_na_pr)

## combine direct and indirect summaries
Both  <-
    mutate(bind_rows(Direct,
                     Indirect),
           qtype = factor(rep(1:2, each = n()/2),
                          labels = c("direct",
                                     "list")))

## difference between list and direct estimates
Diff <- 
    gather(Both,
           varname,
           varval,
           mean:na_pr) %>%
    unite(tempvar,
          qtype,
          varname) %>%
    spread(tempvar,
           varval) %>%
    transmute(year,
              metric,
              faction,
              mean = list_mean - direct_mean,
              se   = sqrt(list_se^2 + direct_se^2),
              p    = round(pnorm(-abs(mean/se)) * 2, 4))

## print out the results:
## "faction" here is actually public religion (low vs high)
## print out the mean and the 95% plus/minus values, along with
## the 95% confidence bounds
## top two rows are for the direct question, bottom two rows
## are for the indirect list question
## top two rows show nearly 100 percent approval regardless
## of level of religiosity
## bottom two rows show 51.7 percent approval for low religiosity,
## and 22.3 percent approval for high religiosity
Both %>%
    filter(!is.na(faction)) %>% 
    mutate(
        pm   = se * 1.96,
        low  = mean - se * 1.96,
        high = mean + se * 1.96
    ) %>%
    select(
        faction,
        mean,
        pm,
        low,
        high
    )

### appendices
#### Appendix B: Sample Composition and Benchmark Comparisons
##### Table 5: demographic benchmarks
###### sources
####### male/female in Lebanon

## UNHCR data source
## https://data2.unhcr.org/en/situations/syria/location/71
## https://data2.unhcr.org/api/population/get/demography?widget_id=71639&geo_id=71&sv_id=4&population_collection=22

## second link is the raw data in JSON format
## it shows that:
##   males:   469150
##   females: 517792
## so the percent male is:
##   469150 / (469150 + 517792) # => .475

####### male/female in Syria

## World Bank data source
## https://datacatalog.worldbank.org/dataset/gender-statistics
## http://databank.worldbank.org/data/reports.aspx?source=gender-statistics

## indicator is:
## Sex ratio at birth (male births per female births)

## available data
## 2002: 1.1
## 2007: 1.1
## 2014: 1.1

####### literacy rates in syria

## World Development Indicators data source

## indicators are
## literacy rate, adult female (% of females ages 15 and above)
## literacy rate, adult male (% of males ages 15 and above)
## literacy rate, adult total (% of people ages 15 and above)

## available data
## 2002: 74.2, 91.0, 82.9
## 2004: 73.6, 87.8, 80.8

## (note the other data source in the AJPS appendix does not
## show later literacy rate data as Erin wrote)

## CIA World Factbook

## indicators are
## literacy, age 15 and over can read and write

## 2015 estimate
## female: 81
## male: 91.7
## total: 86.4

####### educational attainment in Syria

## World Development Indicators data source

## college education
## Educational attainment, at least completed post-secondary, population 25+, female (%) (cumulative)
## Educational attainment, at least completed post-secondary, population 25+, male (%) (cumulative)
## Educational attainment, at least completed post-secondary, population 25+, total (%) (cumulative)

## 2009 estimate
## 10.3, 14.5, 12.4

## HS or better
## Educational attainment, at least completed upper secondary, population 25+, female (%) (cumulative)
## Educational attainment, at least completed upper secondary, population 25+, male (%) (cumulative)
## Educational attainment, at least completed upper secondary, population 25+, total (%) (cumulative)

## 2009 estimate
## 18.9, 24.8, 21.9


## middle school or better
## Educational attainment, at least completed lower secondary, population 25+, female (%) (cumulative)
## Educational attainment, at least completed lower secondary, population 25+, male (%) (cumulative)
## Educational attainment, at least completed lower secondary, population 25+, total (%) (cumulative)

## 2009 estimate
## 29.0, 38.9, 34.1

####### educational attainment in Lebanon

## source
## ASSESSMENT OF THE IMPACT OF SYRIAN REFUGEES IN LEBANON AND THEIR
## EMPLOYMENT PROFILE 2013 International Labour Organization Regional
## Office for the Arab States
## p. 20

####### household size in Syria

## source
## Food security in Syria: Preliminary results based on the
## 2006/07 expenditure survey
## p. 7

####### household size in Lebanon

## source
## ASSESSMENT OF THE IMPACT OF SYRIAN REFUGEES IN LEBANON AND THEIR
## EMPLOYMENT PROFILE 2013 International Labour Organization Regional
## Office for the Arab States
## p. 19

####### ethnic/religious

## CIA World Factbook
## Arab 90.3; Kurdish, Armenian, other 9.7

## Sunni 74
## Alawi, Ismaili, and Shia 13
## Christian 10
## Druze 3

####### poverty

## World Factbook

## population below poverty line (2014 estimate)
## 82.5

###### gather data

## means by year
Syr_main <- 
    summarize(
        select(
            group_by(AC, year),
            male,
            educ_literate,
            educ_hs,
            educ_college,
            syr_fam_size,
            leb_fam_size
        ),
        male = mean(male, na.rm = TRUE),
        literate = mean(educ_literate, na.rm = TRUE),
        hs = mean(educ_hs, na.rm = TRUE),
        college = mean(educ_college, na.rm = TRUE),
        syr = mean(syr_fam_size, na.rm = TRUE),
        leb = mean(leb_fam_size, na.rm = TRUE)
    )

## means by year by sex
Syr_bymale <- 
    summarize(
        select(
            group_by(AC, year, male),
            educ_literate,
            educ_hs,
            educ_college
        ),
        literate = mean(educ_literate, na.rm = TRUE),
        hs = mean(educ_hs, na.rm = TRUE),
        college = mean(educ_college, na.rm = TRUE)
    ) %>%
    gather(varname,
           varval,
           literate:college) %>%
    unite(tempvar,
          male,
          varname) %>%
    spread(tempvar,
           varval) %>%
    rename_at(vars(starts_with("FALSE")),
              funs(str_replace(., "FALSE", "female"))) %>%
    rename_at(vars(starts_with("TRUE")),
              funs(str_replace(., "TRUE", "male")))

## benchmarks in table order
Syr_benchmarks <-
    left_join(Syr_main, Syr_bymale) %>%
    select(
        male,
        female_literate,
        male_literate,
        literate,
        female_hs,
        male_hs,
        hs,
        female_college,
        male_college,
        college,
        syr,
        leb
    ) %>%
    round(digits = 2) %>%
    t() %>%
    as_tibble()
colnames(Syr_benchmarks) <- c("leb15", "leb17")

Benchmarks <-
    tribble(
        ~ stat,               ~ syr07, ~ syr09,  ~ syr15, ~ leb,
        "male",               NA,      .52,      NA,      .48,  
        "literacy, female",   NA,      NA,       .81,     .85,  
        "literacy, male",     NA,      NA,       .92,     .87,  
        "literacy, total",    NA,      NA,       .86,     .86,  
        "HS, female",         NA,      .19,      NA,      .10,  
        "HS, male",           NA,      .25,      NA,      .13,  
        "HS, total",          NA,      .22,      NA,      .12,  
        "college, female",    NA,      .10,      NA,      .03,  
        "college, male",      NA,      .15,      NA,      .04,  
        "college, total",     NA,      .12,      NA,      .04,  
        "family size syr",    5.8,     NA,       NA,      NA,  
        "family size leb",    NA,      NA,       NA,      5
    )                          

Benchmarks <- bind_cols(Benchmarks, Syr_benchmarks)

Benchmarks[,"stat"] <-
    c("Proportion Male",
      "Literacy, Female",
      "Literacy, Male",
      "Literacy, Total",
      "Secondary School Completion, Female",
      "Secondary School Completion, Male",
      "Secondary School Completion, Total",
      "College Completion, Female",
      "College Completion, Male",
      "College Completion, Total",
      "Family Size, Pre-War",
      "Family Size, Post-War")

names(Benchmarks) <-
    c("Statistic",
      "2007",                           # Syria 2007
      "2009",                           # Syria 2009
      "2015",                           # Syria 2015
      "Population",                     # Lebanon population
      "2015 Sample",                    # Lebanon sample 2015
      "2017 Sample")                    # Lebanon sample 2017

###### table

## fancy version
print(xtable(Benchmarks,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 6)),
             digits  = 2,
             label   = "tab:benchmarks",
             caption = "Basic Demographic Benchmarks"),
      table.placement        = "htbp",
      add.to.row = list(pos      = list(0),
                         command = 
                             ## samples, midrules, sumstat names
                             paste0(
                                 paste0("& ",
                                        paste("\\multicolumn{3}{c}{Syria}",
                                              "\\multicolumn{3}{c}{Lebanon}",
                                              sep = " & "),
                                        " \\\\ \n",
                                        collapse = ""),
                                 paste0("\\cmidrule(lr){2-4}",
                                        "\\cmidrule(lr){5-7}",
                                        "\n",
                                        collapse = " "),
                                 paste0("& ", # samples
                                        paste("\\multicolumn{1}{c}{2007}",
                                              "\\multicolumn{1}{c}{2009}",
                                              "\\multicolumn{1}{c}{2015}",
                                              "\\multicolumn{1}{c}{Population}",
                                              "\\multicolumn{1}{c}{2015 Sample}",
                                              "\\multicolumn{1}{c}{2017 Sample}",
                                              sep = " & "),
                                        " \\\\ \n",
                                        collapse = "")
                             )),
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

##### Figure 6: age distributions
###### process data
####### gather data from UN-sponsored Vulnerability Assessment

## UN figures compiled from:
##   UNICEF et al. (2017), Vulnerability Assessment of Syrian
##   Refugees in Lebanon, 2016, produced jointly by UNICEF,
##   UNHCR, and the UN World Food Programme.
##     authors sent us the data file for the age distribution.
##     see <unicef-etal-vulnerability-assessment-age.xlsx>

UN_age <-
    data.frame(bins   = factor(1:15,
                               labels = c("0-4",
                                          "5-9",
                                          "10-14",
                                          "15-19",
                                          "20-24",
                                          "25-29",
                                          "30-34",
                                          "35-39",
                                          "40-44",
                                          "45-49",
                                          "50-54",
                                          "55-59",
                                          "60-64",
                                          "65-69",
                                          "70+")),
               male   = c(20, 18, 13, 8, 5,
                          6,  9,  7,  4, 3,
                          2,  2,  1,  1, 1) / 100,
               female = c(17, 16, 12, 9, 9,
                          10, 8,  6,  4, 3,
                          2,  2,  1,  1, 1) / 100)
UN_age <-
    with(UN_age,
         data.frame(bins = bins,
                    age  = c(male, female),
                    sex  = rep(factor(1:2,
                                      labels = c("Male",
                                                 "Female")),
                               each = nrow(UN_age))))

####### gather survey data

## bin survey data
AC$agebins5 <-
    with(AC,
         cut(age,
             c(-Inf, seq(24, 64, 5), Inf),
             labels = c("20-24",
                        "25-29",
                        "30-34",
                        "35-39",
                        "40-44",
                        "45-49",
                        "50-54",
                        "55-59",
                        "60-64",
                        "65-69")))

## bin survey data to match the UN data
AC$agebins5 <-
    AC$agebins5_alt <-
        with(AC,
             cut(age,
                 c(-Inf, seq(4, 69, 5), Inf),
                 labels = c("0-4",
                            "5-9",
                            "10-14",
                            "15-19",
                            "20-24",
                            "25-29",
                            "30-34",
                            "35-39",
                            "40-44",
                            "45-49",
                            "50-54",
                            "55-59",
                            "60-64",
                            "65-69",
                            "70+")))

Dat15 <- filter(AC, year == 2015)
Dat17 <- filter(AC, year == 2017)

## get the survey data into dataframes (2015)
Age15 <-
    data.frame(bins = levels(Dat15$agebins5),
               age  = c(with(subset(Dat15, female == 0),
                             prop.table(table(agebins5))),
                       with(subset(Dat15, female == 1),
                            prop.table(table(agebins5)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Dat15$agebins5))))
Age15_alt <-
    data.frame(bins = ordered(1:length(levels(Dat15$agebins5_alt)),
                              labels = levels(Dat15$agebins5_alt)),
               age  = c(with(subset(Dat15, female == 0),
                             prop.table(table(agebins5_alt))),
                       with(subset(Dat15, female == 1),
                            prop.table(table(agebins5_alt)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Dat15$agebins5_alt))))

## get the survey data into dataframes (2017)
Age17 <-
    data.frame(bins = levels(Dat17$agebins5),
               age  = c(with(subset(Dat17, female == 0),
                             prop.table(table(agebins5))),
                       with(subset(Dat17, female == 1),
                            prop.table(table(agebins5)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Dat17$agebins5))))
Age17_alt <-
    data.frame(bins = ordered(1:length(levels(Dat17$agebins5_alt)),
                              labels = levels(Dat17$agebins5_alt)),
               age  = c(with(subset(Dat17, female == 0),
                             prop.table(table(agebins5_alt))),
                       with(subset(Dat17, female == 1),
                            prop.table(table(agebins5_alt)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Dat17$agebins5_alt))))

####### put UN and survey data together (adults only)

## UN data adults only
UN_adult <-
    subset(UN_age,
           !(bins %in% c("0-4",
                         "5-9",
                         "10-14",
                         "15-19")))
UN_adult$age <- 
    c(with(subset(UN_adult, sex == "Male"),
           age / sum(age)),
      with(subset(UN_adult, sex == "Female"),
           age / sum(age)))

## survey data adults only
Survey15_adult <- 
    subset(Age15_alt,
           !(bins %in% c("0-4",
                         "5-9",
                         "10-14",
                         "15-19")))
Survey17_adult <- 
    subset(Age17_alt,
           !(bins %in% c("0-4",
                         "5-9",
                         "10-14",
                         "15-19")))


## combine UN and survey data
Adult_combo <- rbind(UN_adult, Survey15_adult, Survey17_adult)
Adult_combo$src <-
    rep(factor(1:3,
               labels = c("UN Vulnerability Assessment",
                          "2015 Survey",
                          "2017 Survey")),
        each = nrow(Adult_combo) / 3)
Adult_combo$age <-
    with(Adult_combo,
         ifelse(sex == "Male", age, -age))

Adult_combo$bins <- factor(Adult_combo$bins)

###### Kolmogorov-Smirnov test of equality of distributions across data sources

## equality of male subsample across survey and UN data
ks.test(filter(Adult_combo,
               sex == "Male",
               src == "UN Vulnerability Assessment")$age,
        filter(Adult_combo,
               sex == "Male",
               src == "2015 Survey")$age)
ks.test(filter(Adult_combo,
               sex == "Male",
               src == "UN Vulnerability Assessment")$age,
        filter(Adult_combo,
               sex == "Male",
               src == "2017 Survey")$age)
ks.test(filter(Adult_combo,
               sex == "Male",
               src == "2015 Survey")$age,
        filter(Adult_combo,
               sex == "Male",
               src == "2017 Survey")$age)

## equality of female subsample across survey and UN data
ks.test(filter(Adult_combo,
               sex == "Female",
               src == "UN Vulnerability Assessment")$age,
        filter(Adult_combo,
               sex == "Female",
               src == "2015 Survey")$age)
ks.test(filter(Adult_combo,
               sex == "Female",
               src == "UN Vulnerability Assessment")$age,
        filter(Adult_combo,
               sex == "Female",
               src == "2017 Survey")$age)
ks.test(filter(Adult_combo,
               sex == "Female",
               src == "2015 Survey")$age,
        filter(Adult_combo,
               sex == "Female",
               src == "2017 Survey")$age)

###### plot UN and survey

## plot it
plot_age_all3 <-
    ggplot(data = Adult_combo,
           aes(x    = bins,
               y    = age,
               fill = sex)) +
    geom_hline(yintercept = 0) +
    geom_bar(stat = "identity",
             data = subset(Adult_combo, sex == "Male"),
             alpha = 0.50) +
    geom_bar(stat = "identity",
             data = subset(Adult_combo, sex == "Female"),
             alpha = 0.50) +
    annotate(geom  = "rect",
             xmin  = nlevels(Adult_combo$bins) - .75,
             xmax  = nlevels(Adult_combo$bins) + .75,
             ymin  = 0.125 * c(-1, 1) - 0.05,
             ymax  = 0.125 * c(-1, 1) + 0.05,
             color = "grey95",
             fill  = "white") +
    annotate(geom     = "text",
             x        = nlevels(Adult_combo$bins),
             y        = -0.125,
             size     = 2,
             fontface = "italic",
             color    = "grey50",
             label    = "Female") +
    annotate(geom     = "text",
             x        = nlevels(Adult_combo$bins),
             y        = 0.125,
             size     = 2,
             fontface = "italic",
             color    = "grey50",
             label    = "Male") +
    scale_fill_manual(values = c("darkred", "darkblue")) +
    scale_y_continuous(limits = c(-0.25, 0.25),
                       breaks = seq(-0.3, 0.3, 0.1),
                       labels = abs(round(seq(-0.3, 0.3, 0.1), 1))) +
    labs(x    = "Age",
         y    = "Proportion") +
    coord_flip() +
    facet_wrap(~ src) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          strip.text       = element_text(size = 8),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_age_all3_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_age_all3),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(1.5, "in")
                         ))

## save it
ggsave(plot     = plot_age_all3_egg,
       filename = "benchmark-age.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 3.00)

##### Figure 7: province of origin
###### gather data

## UNHCR figures compiled from:
##   UNHCR Lebanon - Beirut Country Office. Syria Refugee Response
##   Lebanon: Places of Origin of Syrian Refugees Registered in
##   Lebanon. 31 March 2015.

## figures on province population in Syria come from:
##   OCHA (2014), "Syrian Arab Republic Governorates Profile (June 2014)"
##   https://reliefweb.int/sites/reliefweb.int/files/resources/Syria%20governorate%20profiles%206%20August%202014.pdf

Province <-
    data.frame(province = levels(AC$province_syr),
               syria_n  = c("Daraa"        = 1027000,
                            "Quneitra"     = 90000,
                            "Al-Suwayda"   = 370000,
                            "Dimashq"      = 1754000,
                            "Rif Dimashq"  = 2835900,
                            "Homs"         = 1803000,
                            "Hama"         = 1628000,
                            "Tartus"       = 797000,
                            "Latakia"      = 1008000,
                            "Idlib"        = 1501000,
                            "Halab"        = 4867991,
                            "Al-Raqqa"     = 944000,
                            "Deir al-Zour" = 1239005,
                            "Al-Hasaka"    = 1512000),
               leb15_n  = as.vector(table(filter(AC, year == 2015)$province_syr)),
               leb17_n  = as.vector(table(filter(AC, year == 2017)$province_syr)),
               ## survey_n = as.vector(table(AC$province_syr)),
               unhcr_n  = c("Daraa"        = 23386 + 30613 + 28354,
                            "Quneitra"     = 10512 + 31,
                            "Al-Suwayda"   = 251 + 666 + 124,
                            "Dimashq"      = 55105,
                            "Rif Dimashq"  = 19175 + 62755 + 14502 + 25305 + 12524 + 2679 + 11440 + 2326,
                            "Homs"         = 1089 + 3072 + 151189 + 15150 + 55705 + 19940,
                            "Hama"         = 16048 + 50585 + 7681 + 11229 + 956,
                            "Tartus"       = 1749 + 25 + 49 + 1181 + 345,
                            "Latakia"      = 3464 + 919 + 40 + 520,
                            "Idlib"        = 4412 + 67666 +10462 + 24061 + 45247,
                            "Halab"        = 4618 + 23882 + 26423 + 28898 + 11140 + 10608 + 133659 + 3094,
                            "Al-Raqqa"     = 5678 + 43411 + 15850,
                            "Deir al-Zour" = 17169 + 3395 + 5094,
                            "Al-Hasaka"    = 2711 + 15233 + 12351 + 618))

## convert to proportions
Province$syria_prop <-
    with(Province,
         syria_n / sum(syria_n))
Province$leb15_prop <-
    with(Province,
         leb15_n / sum(leb15_n))
Province$leb17_prop <-
    with(Province,
         leb17_n / sum(leb17_n))
Province$unhcr_prop <-
    with(Province,
         unhcr_n / sum(unhcr_n))
## difference between survey proportions and UNHCR proportions
Province$leb15_vs_unhcr <-
    with(Province,
         leb15_prop - unhcr_prop)
Province$leb17_vs_unhcr <-
    with(Province,
         leb17_prop - unhcr_prop)
## fancy province labels: province name + proportion of Syrian pop
levels(Province$province) <-
    with(Province,
         sort(paste(as.character(province),
                    paste0("(",
                           substr(format(round(Province$syria_prop, 2)), 2, 10),
                           ")"))))
## reorder province levels by survey proportion
Province$province <-
    with(Province,
         reorder(province, -unhcr_prop))

## get a quick summary of the proportions by sample and UNHCR
Province_3val <-
    select(Province,
           province,
           leb15_prop,
           leb17_prop,
           unhcr_prop) %>%
    gather(source,
           value,
           leb15_prop:unhcr_prop)

Province_3val$source <-
    with(Province_3val,
         factor(source,
                labels = c("2015 Sample",
                           "2017 Sample",
                           "UNHCR")))

###### Kolmogorov-Smirnov  test of equality of distributions across data sources

with(Province,
     ks.test(leb15_prop, unhcr_prop))

with(Province,
     ks.test(leb17_prop, unhcr_prop))

with(Province,
     ks.test(leb15_prop, leb17_prop))

###### plot

## plot it
plot_province <-
    ggplot(data = Province_3val,
           aes(x    = province,
               y    = value,
               fill = source)) +
    geom_bar(position = "dodge",
             stat     = "identity") +
    coord_flip() +
    labs(x    = "Province",
         y    = "Proportion") +
    scale_fill_manual(values = c("grey75", "grey65", "dodgerblue4")) +
    guides(fill=guide_legend(title = "Data Source")) +
    theme_bw() +
    theme(legend.position  = "bottom",
          strip.background = element_blank(),
          strip.text       = element_text(size = 8),
          panel.border     = element_rect(color = "black"))
plot_province

## save it
ggsave(plot     = plot_province,
       filename = "benchmark-province.pdf",
       family   = "Palatino",
       width    = 6,
       height   = 6)

##### Figure 8: province/region in Lebanon
###### gather data

## UNHCR figures from:
##  Situation Syria Regional Refugee Response
##  https://data2.unhcr.org/en/situations/syria/location/71
##    data from wikipedia page:
##      https://en.wikipedia.org/wiki/Syrians_in_Lebanon
##      retrieved 2018-12-19
##    I'm using the wikipedia data because the UNHCR figures
##    are updated frequently and I don't want to go dig around
##    on their website to find 2017-era data.  The figures don't
##    change dramatically over time, but late-2018 is closer to
##    the sample than is late-2020 when I'm making this figure.

Lebprov <- 
    tibble(province =
               c("Mount Lebanon",
                 "North",
                 "South",
                 "Bekaa"),
           leb = 
               prop.table(table(AC$province_leb_4))[c("Mount Lebanon",
                                                      "North",
                                                      "South",
                                                      "Bekaa")],
           ab =
               prop.table(
                   table(
                       factor(
                           filter(ABSY, lebanon)$province_4)
                   )
               )[c("Mount Lebanon",
                   "North",
                   "South",
                   "Bekaa")],
           unhcr =
               c(.261,
                 .261,
                 .120,
                 .358)
           )

## reorder
Lebprov$province <-
    with(Lebprov,
         reorder(province, -unhcr))

## get a quick summary of the proportions by sample and UNHCR
Lebprov_3val <-
    select(Lebprov,
           province,
           leb,
           ab,
           unhcr) %>%
    gather(source,
           value,
           leb:unhcr)
Lebprov_3val$source <-
    with(Lebprov_3val,
         factor(source,
                labels = c("Original Sample",
                           "Arab Barometer",
                           "UNHCR")))

###### Kolmogorov-Smirnov test of equality of distributions across data sources

with(Lebprov,
     ks.test(leb, ab))

with(Lebprov,
     ks.test(leb, unhcr))

with(Lebprov,
     ks.test(ab, unhcr))

###### plot

## plot it
plot_province <-
    ggplot(data = Lebprov_3val,
           aes(x    = province,
               y    = value,
               fill = source)) +
    geom_bar(position = "dodge",
             stat     = "identity") +
    coord_flip() +
    labs(x    = "Province",
         y    = "Proportion") +
    scale_fill_manual(values = c("grey75", "grey65", "dodgerblue4")) +
    guides(fill=guide_legend(title = "Data Source")) +
    theme_bw() +
    theme(legend.position  = "bottom",
          strip.background = element_blank(),
          strip.text       = element_text(size = 8),
          panel.border     = element_rect(color = "black"))

## save it
ggsave(plot     = plot_province,
       filename = "benchmark-province-leb.pdf",
       family   = "Palatino",
       width    = 6,
       height   = 6)

#### Appendix D: Rank Order Question
##### Table 6: factional rank order profiles
###### helper function

## helper function: recode NAs as "." to make tables easier to read
clean_factions <- function(x) {
    x <- factor(ifelse(x < 0 | is.na(x), 4, x))
    x <- recode(x, "4" = ".")
}

###### process data
####### 2015

## clean up the raw faction responses
Fac15 <-
    Syr15 %>%
    transmute(
        fsa       = clean_factions(Q28a),
        govt      = clean_factions(Q28b),
        islamistd = clean_factions(Q28c),
        islamistf = clean_factions(Q28d),
        kurds     = clean_factions(Q28e),
        hizballah = clean_factions(Q28f)
    )

## get each unique response profile and count them up
## https://stackoverflow.com/questions/51272510/how-to-count-unique-rows-in-a-data-frame
Profiles15 <- 
    Fac15 %>%
    select(
        `Syrian government` = govt,
        `Hizballah`         = hizballah,
        `Kurdish groups`    = kurds,
        `Free Syrian Army`  = fsa,
        `Syrian Islamists`  = islamistd,
        `Foreign Islamists` = islamistf
    ) %>%
    group_by_all %>%
    count

####### 2017

## clean up the raw faction responses
Fac17 <-
    Syr17 %>%
    transmute(
        fsa       = clean_factions(Q24a),
        govt      = clean_factions(Q24b),
        islamistd = clean_factions(Q24c),
        islamistf = clean_factions(Q24d),
        kurds     = clean_factions(Q24e),
        hizballah = clean_factions(Q24f)
    )

## get each unique response profile and count them up
## https://stackoverflow.com/questions/51272510/how-to-count-unique-rows-in-a-data-frame
Profiles17 <- 
    Fac17 %>%
    select(
        `Syrian government` = govt,
        `Hizballah`         = hizballah,
        `Kurdish groups`    = kurds,
        `Free Syrian Army`  = fsa,
        `Syrian Islamists`  = islamistd,
        `Foreign Islamists` = islamistf
    ) %>%
    group_by_all %>%
    count

###### print

## 2015
print(xtable(Profiles15,
             align   = c("l ",
                         rep("C{6em} ", 6),
                         "r"),
             caption = "Factional Rank Order Profiles, 2015",
             label   = "tab:profiles:2015"),
      table.placement = "htbp",
      scalebox        = 0.75,
      booktabs        = TRUE)

## 2017
print(xtable(Profiles17,
             align   = c("l ",
                         rep("C{6em} ", 6),
                         "r"),
             caption = "Factional Rank Order Profiles, 2017",
             label   = "tab:profiles:2017"),
      table.placement = "htbp",
      scalebox        = 0.75,
      booktabs        = TRUE)

##### Table 7: faction by Lebanese region
###### gather data

Dat0 <-
    transmute(
        AC,
        fac3 =
            factor(ifelse(fac_4cat == "none", 3,
                   ifelse(fac_4cat == "govt", 1,
                   ifelse(fac_4cat == "oppn", 2,
                          NA))),
                   labels = c("Government",
                              "Rebels",
                              "Unaligned")),
        fac3 =
            factor(fac3,
                   labels =
                       paste0(
                           levels(fac3),
                           " (",
                           round(prop.table(table(fac3)) * 100, 1) ,
                           ")"
                       )
                   ),
        prov4 = province_leb_4,
        prov4 =
            factor(prov4,
                   labels =
                       paste0(
                           levels(prov4),
                           " (",
                           round(prop.table(table(prov4)) * 100, 1) ,
                           ")"
                       )
                   )
    )

###### print

provtab <- 
    with(
        Dat0,
        round(prop.table(table(fac3, prov4), margin = 2) * 100, 1)
    )

## print it out
print(xtable(provtab,
             align   = c ("l ",
                          ## "l ", # extra for (removed) row numbers
                          rep("D{.}{.}{2.2} ", 4)),
             digits  = 1,
             caption = paste("Factional Support in Each Lebanese Region.",
                             "Marginal distributions of faction and region",
                             "reported in parentheses.",
                             "Columns may not add to 100 due to rounding."
                             ),
             label   = "tab:lebprovince"
             ),
      table.placement = "htbp",
      ## ugh: getting a pretty table header is a mess
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # column headers
                                         paste0("\\multicolumn{1}{c}{",
                                                colnames(provtab),
                                                "}",
                                                collapse = " & "),
                                         " \\\\ \n")),
      include.rownames = TRUE,
      include.colnames = FALSE,
      booktabs         = TRUE
      )

##### Table 8: rank order non-alignment

## model no faction
mod_rank_none <-
    glm(fac_none
        ~ young
        * survey17
        + male
        * survey17
        + minority
        * survey17
        + log(syr_room_density)
        * survey17
        + log(leb_room_density)
        * survey17
        + educ_hs
        * survey17
        + interest_bi
        * survey17
        + quran_weekly
        * survey17
        + relpol_bi
        * survey17
       ,
        family = binomial,
        data   = AC)

tex_mods <- list(mod_rank_none)

tex_cols <- c("Unaligned")

## variable names (original order)
tex_vars <-
    c("(Intercept)",
      "Young",
      "2017",
      "Male",
      "Minority",
      "Crowded in Syria",
      "Crowded in Lebanon",
      "Educated",
      "Interested in Politics",
      "Private Religion",
      "Public Religion",
      "2017 $\\times$ Older",
      "2017 $\\times$ Female",
      "2017 $\\times$ Minority",
      "2017 $\\times$ Crowded in Syria",
      "2017 $\\times$ Crowded in Lebanon",
      "2017 $\\times$ Educated",
      "2017 $\\times$ Interested in Politics",
      "2017 $\\times$ Private Religion",
      "2017 $\\times$ Public Religion")

## reorder coefficients
tex_cord <- c(1:2, 4:11, 3, 12:20)

## texreg it to screen
texreg(l                  = tex_mods,
       custom.coef.names  = tex_vars,
       custom.model.names = tex_cols,
       reorder.coef       = tex_cord,
       single.row         = TRUE,
       leading.zero       = TRUE,
       booktabs           = TRUE,
       dcolumn            = TRUE,
       stars              = c(.01, .05))

#### Appendix E: Balance Checks
##### Table 9: balance checks
###### gather data

## container
balance <- list()

## gather data
Dat0 <-
    transmute(
        AC,
        year       = year,
        group      = ac_group,
        age        = age,
        female     = female,
        minority   = minority,
        educ       = educ_hs,
        lnsrd      = log(syr_room_density),
        lnlrd      = log(leb_room_density),
        know       = know_bi,
        interest   = interest_bi,
        understand = understand_bi,
        engaged    = engaged_bi,
        pray       = pray_daily,
        quran      = quran_daily,
        rpol       = rel_pol,
        recon      = rel_econ,
        rcrime     = rel_crime,
        rfam       = rel_family,
        relpol     = relpol_bi,
        fac        = fac_all
    )

## pretty variable names
varlabs  <-
    c("Age",
      "Female",
      "Minority",
      "Education",
      "Room Density in Syria",
      "Room Density in Lebanon",
      "Knowledge of Politics",
      "Interest in Politics",
      "Understand Politics",
      "Engaged in Politics",
      "Prayer",
      "Quran/Bible",
      "Religion in Politics",
      "Religion in Economy",
      "Religion in Crime",
      "Religion in Family",
      "Religion Index",
      "Faction")

## pretty summary names
sumlabs <-
    c("Minimum",
      "1st Quarter",
      "Median",
      "3rd Quarter",
      "Maximum")

## loop through the samples (ugh)
for (y in c(2015, 2017)) {
    Dat00    <- filter(Dat0, year == y)
    varnames <- names(select(Dat0, -year, -group))
    bal <- 
        matrix(data     = NA,
               nrow     = length(varnames),
               ncol     = 2,
               dimnames = list(varnames,
                               c("stat", "p")))
    for (covariate in varnames) {
        ## categorical if < 10 values, else continuous
        if (with(Dat00,
                 length(unique(eval(parse(text = covariate)))) < 10)) {
            ## categorical: use chi-square test
            test <-
                with(Dat00,
                     chisq.test(table(group,
                                      eval(parse(text = covariate)))))
        } else {
            ## continuous: use F-test
            test <- 
                with(Dat00,
                     summary(aov(eval(parse(text = covariate))
                                 ~ group)))
            test$statistic <- test[[1]][[1, "F value"]]
            test$p.value   <- test[[1]][[1, "Pr(>F)"]]
        }
        bal[covariate,] <-
            c(test$statistic,
              test$p.value)
    }
    bal <- cbind(bal, "fdr" = p.adjust(bal[,"p"], method = "fdr"))
    balance[[paste0("year", y)]] <- bal
}

## get the balance p-values for 2015
Bal15 <-
    with(balance,
         tibble(var = c(varlabs, sumlabs),
                p   = c(year2015[,"p"],
                        quantile(year2015[,"p"],
                                 c(0, .25, .50, .75, 1))),
                fdr = c(year2015[,"fdr"],
                        quantile(year2015[,"fdr"],
                                 c(0, .25, .50, .75, 1)))))

## get the balance p-values for 2017
Bal17 <-
    with(balance,
         tibble(var = c(varlabs, sumlabs),
                p   = c(year2017[,"p"],
                        quantile(year2017[,"p"],
                                 c(0, .25, .50, .75, 1))),
                fdr = c(year2017[,"fdr"],
                        quantile(year2017[,"fdr"],
                                 c(0, .25, .50, .75, 1)))))

## combine the balance values
Bal <-
    mutate(bind_rows(Bal15,
                     Bal17),
           sample = rep(factor(1:2,
                               labels = c("2015",
                                          "2017")),
                        each = n() / 2)) %>%
    gather(varname,
           varval,
           p:fdr) %>%
    unite(tempvar,
          varname,
          sample) %>%
    spread(tempvar,
           varval) %>%
    select(var,
           p_2015,
           fdr_2015,
           p_2017,
           fdr_2017) %>%
    ## select alphabetizes var, so we need to retrieve the right order
    arrange(match(var,
                  c(varlabs,
                    sumlabs)))

###### print

## print it out
print(xtable(Bal,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 4)),
             caption = "Balance Checks.  Raw and False Discovery Rate-corrected $p$-values.",
             label   = "tab:balance"),
      table.placement  = "htbp",
      ## ugh: getting a pretty table header is a mess
      add.to.row       = list(pos     = list(0, 0, 0, nrow(Bal) - 5),
                              command =
                                  ## samples, midrules, sumstat names
                                  c(paste0("& ", # samples
                                           paste("\\multicolumn{2}{c}{2015}",
                                                 "\\multicolumn{2}{c}{2017}",
                                                 sep = " & "),
                                           " \\\\ \n",
                                           collapse = ""),
                                    paste(paste("\\cmidrule(lr){2-3}", # midrules
                                                "\\cmidrule(lr){4-5}",
                                                collapse = " & "),
                                          "\n"),
                                    paste0("& ", # counts
                                           paste("\\multicolumn{1}{c}{Raw}",
                                                 "\\multicolumn{1}{c}{FDR}",
                                                 "\\multicolumn{1}{c}{Raw}",
                                                 "\\multicolumn{1}{c}{FDR}",
                                                 sep = " & "),
                                           " \\\\ \n",
                                           collapse = ""),
                                    "\\midrule \n")),
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

#### Appendix F: Non-Sensitive Items and Raw List Counts
##### Table 10: support for non-sensitive items
###### gather data

## gather data
Dat0 <-
    transmute(
        AC,
        sample =
            factor(ifelse(year == 2015, 1, 2),
                   labels = c("y2015",
                              "y2017")),
        control    = ac_group_control,
        parl       = ac_parl_direct,
        noparty    = ac_no_parties_direct,
        mil        = ac_mil_direct,
        pres       = ac_pres_direct,
        multiparty = ac_multiparty_direct,
        mandate    = ac_mandate_direct
    )

## means of each item by sample as *percentage*
offitems_1517 <-
    summarize_all(
        group_by(filter(Dat0, control),
                 sample),
        funs(mean = mean(., na.rm = TRUE) * 100)
    )

## means of each item pooled across samples as *percentage*
offitems_pool <-
    summarize_all(
        group_by(
            mutate(
                filter(Dat0, control),
                sample = factor(3, labels = "pooled")
            ),
            sample),
        funs(mean = mean(., na.rm = TRUE) * 100)
    )

## combine the summaries
offitems <- 
    select(bind_rows(offitems_pool,
                     offitems_1517),
           -control_mean) %>%
    gather(key   = varname,
           value = varval,
           2:ncol(.)) %>%
    spread(key   = names(.)[1],
           value = varval) %>%
    arrange(match(varname,
                  c("parl_mean",
                    "noparty_mean",
                    "mil_mean",
                    "pres_mean",
                    "multiparty_mean",
                    "mandate_mean")))

offitems$varname <-
    c("Parliamentary democracy",
      "No-party system",
      "Military government",
      "Presidential democracy",
      "Multiparty system",
      "French Mandate")

###### print

## print it out
print(xtable(offitems,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 3)),
             digits  = 1,
             caption = "Percent Support for Non-Sensitive Options",
             label   = "tab:offitems"),
      table.placement  = "htbp",
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # samples
                                         paste("\\multicolumn{1}{c}{Pooled}",
                                               "\\multicolumn{1}{c}{2015}",
                                               "\\multicolumn{1}{c}{2017}",
                                               sep = " & "),
                                         " \\\\ \n",
                                         collapse = "")),
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

##### Table 11: item counts, Assad lists
###### gather data

## gather data
Dat0 <-
    transmute(
        AC,
        sample = year,
        treat  = factor(ifelse(ac_group_control, 1, 2),
                        labels = c("Control",
                                   "Treatment")),
        count  = ac_assad_list
    )

## 2015 in *percentage* terms
rawprops15 <- 
    as_tibble(with(filter(Dat0, sample == 2015),
                   prop.table(table(treat,
                                    count),
                              margin = 1) * 100)) %>%
    spread(count, n)
rawprops15$`4` <- c(NA, 0)              # NA for control 4 count

## 2017 in *percentage* terms
rawprops17 <- 
    as_tibble(with(filter(Dat0, sample == 2017),
                   prop.table(table(treat,
                                    count),
                              margin = 1) * 100)) %>%
    spread(count, n)                    # NA for control 4 count
rawprops17$`4`[1] <- NA

###### print

## print it out: 2015
print(xtable(rawprops15,
             digits  = 1,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 5)),
             caption = "Raw Count Percentage, Assad 2015",
             label   = "tab:rawcount:assad:2015"),
      table.placement  = "htbp",
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # counts
                                         paste("\\multicolumn{1}{c}{0}",
                                               "\\multicolumn{1}{c}{1}",
                                               "\\multicolumn{1}{c}{2}",
                                               "\\multicolumn{1}{c}{3}",
                                               "\\multicolumn{1}{c}{4}",
                                               sep = " & "),
                                         " \\\\ \n",
                                         collapse = "")),
      NA.string              = "",
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

## print it out: 2017
print(xtable(rawprops17,
             digits  = 1,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 5)),
             caption = "Raw Count Percentage, Assad 2017",
             label   = "tab:rawcount:assad:2017"),
      table.placement  = "htbp",
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # counts
                                         paste("\\multicolumn{1}{c}{0}",
                                               "\\multicolumn{1}{c}{1}",
                                               "\\multicolumn{1}{c}{2}",
                                               "\\multicolumn{1}{c}{3}",
                                               "\\multicolumn{1}{c}{4}",
                                               sep = " & "),
                                         " \\\\ \n",
                                         collapse = "")),
      NA.string              = "",
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

##### Table 12: item counts, jihadi lists
###### gather data

## gather data
Dat0 <-
    transmute(
        AC,
        sample = year,
        treat  = factor(ifelse(ac_group_control, 1, 2),
                        labels = c("Control",
                                   "Treatment")),
        count  = ac_caliphate_list
    )

## 2015 in *percentage* terms
rawprops15 <- 
    as_tibble(with(filter(Dat0, sample == 2015),
                   prop.table(table(treat,
                                    count),
                              margin = 1) * 100)) %>%
    spread(count, n)
rawprops15$`3` <- c(0,  0)              # no one counts 3
rawprops15$`4` <- c(NA, 0)              # NA for control 4 count

## 2017 in *percentage* terms
rawprops17 <- 
    as_tibble(with(filter(Dat0, sample == 2017),
                   prop.table(table(treat,
                                    count),
                              margin = 1) * 100)) %>%
    spread(count, n)                    # NA for control 4 count
rawprops17$`4` <- c(NA, 0)              # NA for control 4 count

###### print

## print it out: 2015
print(xtable(rawprops15,
             digits  = 1,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 5)),
             caption = "Raw Count Percentage, Jihadi 2015",
             label   = "tab:rawcount:jihad:2015"),
      table.placement  = "htbp",
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # counts
                                         paste("\\multicolumn{1}{c}{0}",
                                               "\\multicolumn{1}{c}{1}",
                                               "\\multicolumn{1}{c}{2}",
                                               "\\multicolumn{1}{c}{3}",
                                               "\\multicolumn{1}{c}{4}",
                                               sep = " & "),
                                         " \\\\ \n",
                                         collapse = "")),
      NA.string              = "",
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

## print it out: 2017
print(xtable(rawprops17,
             digits  = 1,
             align   = c("l ",
                         "l ", # extra for (removed) row numbers
                         rep("D{.}{.}{2.2} ", 5)),
             caption = "Raw Count Percentage, Jihadi 2017",
             label   = "tab:rawcount:jihad:2017"),
      table.placement  = "htbp",
      add.to.row       = list(pos     = list(0),
                              command =
                                  paste0("& ", # counts
                                         paste("\\multicolumn{1}{c}{0}",
                                               "\\multicolumn{1}{c}{1}",
                                               "\\multicolumn{1}{c}{2}",
                                               "\\multicolumn{1}{c}{3}",
                                               "\\multicolumn{1}{c}{4}",
                                               sep = " & "),
                                         " \\\\ \n",
                                         collapse = "")),
      NA.string              = "",
      include.rownames       = FALSE,
      include.colnames       = FALSE,
      sanitize.text.function = function(x){x},
      booktabs               = TRUE)

#### Appendix G: List Experiment Results Split by Metric
##### Figure 9: support by metric and year
###### Assad

## gather list experiment data
Dat0 <-
    transmute(
        Assad$base$main$full,
        mean,
        se,
        year = year,
        metric = factor(metric,
                        levels = levels(metric)[2:1],
                        labels = c("Okay",
                                   "Good")),
        qtype = factor(qtype,
                       labels = c("Direct",
                                  "List"))
    )

## plot it
plot_assad_metricyear <-
    ggplot(Dat0,
           aes(x    = metric,
               y    = mean)) +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.15, 0.65),
                       labels       = seq(0.0, 1.0, .25) * 100,
                       breaks       = seq(0.0, 1.0, .25),
                       minor_breaks = seq(-0.125, 0.825, .25)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_grid(year ~ qtype) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_assad_metricyear_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_assad_metricyear),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_assad_metricyear_egg,
       filename = "assad-splitmetrics-1517.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 2.00)

###### Jihad

## gather list experiment data
Dat0 <-
    transmute(
        Jihad$base$main$full,
        mean,
        se,
        year = year,
        metric = factor(metric,
                        levels = levels(metric)[2:1],
                        labels = c("Okay",
                                   "Good")),
        qtype = factor(qtype,
                       labels = c("Direct",
                                  "List"))
    )

## plot it
plot_jihad_metricyear <-
    ggplot(Dat0,
           aes(x    = metric,
               y    = mean)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.30, 0.30),
                       labels       = seq(-1.0, 1.0, .20) * 100,
                       breaks       = seq(-1.0, 1.0, .20),
                       minor_breaks = seq(-1.0, 1.0, .10)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_grid(year ~ qtype) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_jihad_metricyear_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_jihad_metricyear),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_jihad_metricyear_egg,
       filename = "jihadi-splitmetrics-1517.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 2.00)

#### Appendix H: List Experiment Results Over Time
##### Figure 10: Assad by factional preference and survey year

## gather list experiment data
Dat0 <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct")),
        ## highlight facets that have direct/indirect differences
        high    = ifelse(faction == "Rebels",
                         FALSE,
                         TRUE)
    )

## plot it
plot_assad_poolmetrics <-
    ggplot(Dat0,
           aes(x    = qtype,
               y    = mean,
               fill = high)) +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## highlight facets
    geom_rect(xmin  = -Inf,
              xmax  =  Inf,
              ymin  = -Inf,
              ymax  =  Inf,
              alpha = 0.035) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.15, 1.15),
                       labels       = seq(0.0, 1.0, .50) * 100,
                       breaks       = seq(0.0, 1.0, .50),
                       minor_breaks = seq(0.25, 0.75, .25)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_grid(year ~ faction) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_assad_poolmetrics_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_assad_poolmetrics),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_assad_poolmetrics_egg,
       filename = "assad-poolmetrics-1517.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 2.00)

##### Figure 11: Jihadis by factional preference and survey year

## gather list experiment data
Dat0 <-
    transmute(
        filter(Jihad$pool$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct")),
        ## highlight estimates with negative probabilities
        neg     = ifelse(mean < 0 &
                         pnorm(-abs(mean/se)) * 2 < .15,
                         TRUE, FALSE),
        ## highlight facets that have direct/indirect differences
        high    = ifelse(faction == "Rebels",
                         FALSE,
                         TRUE)
    )

## plot it
plot_jihad_poolmetrics <-
    ggplot(Dat0,
           aes(x     = qtype,
               y     = mean,
               color = neg,
               fill  = high)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## highlight facets
    geom_rect(xmin  = -Inf,
              xmax  =  Inf,
              ymin  = -Inf,
              ymax  =  Inf,
              alpha = 0.035) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.30, 0.30),
                       labels       = seq(-1.0, 1.0, .20) * 100,
                       breaks       = seq(-1.0, 1.0, .20),
                       minor_breaks = seq(-1.0, 1.0, .10)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    scale_color_manual(values = c("black", "red4")) +
    coord_flip() +
    facet_grid(year ~ faction) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_jihad_poolmetrics_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_jihad_poolmetrics),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_jihad_poolmetrics_egg,
       filename = "jihadi-poolmetrics-1517.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 2.00)

##### Figure 12: chance in support over time
###### gather Assad list data

## gather list experiment data
Dat_list <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               qtype   == "list",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )

Dat_direct <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               qtype   == "direct",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )


## ugh, break this up
L_Govt <- subset(Dat_list, faction == "Government")
L_Rebs <- subset(Dat_list, faction == "Rebels")
L_None <- subset(Dat_list, faction == "Unaligned")

D_Govt <- subset(Dat_direct, faction == "Government")
D_Rebs <- subset(Dat_direct, faction == "Rebels")
D_None <- subset(Dat_direct, faction == "Unaligned")

L_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(L_Govt$mean * c(-1, 1)),
                     sum(L_Rebs$mean * c(-1, 1)),
                     sum(L_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(L_Govt$se^2)),
                     sqrt(sum(L_Rebs$se^2)),
                     sqrt(sum(L_None$se^2))))

D_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(D_Govt$mean * c(-1, 1)),
                     sum(D_Rebs$mean * c(-1, 1)),
                     sum(D_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(D_Govt$se^2)),
                     sqrt(sum(D_Rebs$se^2)),
                     sqrt(sum(D_None$se^2))))

Assad_Diff <-
    mutate(
        bind_rows(D_Diff,
                  L_Diff),
        qtype =
            factor(rep(2:1, each = 3),
                   labels = c("List",
                              "Direct")),
        high_point =
            ifelse((faction == "Government" & qtype == "List") |
                   (faction == "Unaligned"  & qtype == "Direct"),
                   TRUE,
                   FALSE),
        high =
            ifelse(faction == "Rebels",
                   FALSE,
                   TRUE)
    )

###### gather Jihadi list data

## gather list experiment data
Dat_list <-
    transmute(
        filter(Jihad$pool$main$sub0, # drop "other"
               qtype   == "list",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )

Dat_direct <-
    transmute(
        filter(Jihad$pool$main$sub0, # drop "other"
               qtype   == "direct",
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned"))
    )


## ugh, break this up
L_Govt <- subset(Dat_list, faction == "Government")
L_Rebs <- subset(Dat_list, faction == "Rebels")
L_None <- subset(Dat_list, faction == "Unaligned")

D_Govt <- subset(Dat_direct, faction == "Government")
D_Rebs <- subset(Dat_direct, faction == "Rebels")
D_None <- subset(Dat_direct, faction == "Unaligned")

L_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(L_Govt$mean * c(-1, 1)),
                     sum(L_Rebs$mean * c(-1, 1)),
                     sum(L_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(L_Govt$se^2)),
                     sqrt(sum(L_Rebs$se^2)),
                     sqrt(sum(L_None$se^2))))

D_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(D_Govt$mean * c(-1, 1)),
                     sum(D_Rebs$mean * c(-1, 1)),
                     sum(D_None$mean * c(-1, 1))),
               se =
                   c(sqrt(sum(D_Govt$se^2)),
                     sqrt(sum(D_Rebs$se^2)),
                     sqrt(sum(D_None$se^2))))

Jihad_Diff <-
    mutate(
        bind_rows(D_Diff,
                  L_Diff),
        qtype =
            factor(rep(2:1, each = 3),
                   labels = c("List",
                              "Direct")),
        high_point =
            ifelse((faction == "Government" & qtype == "List") |
                   (faction == "Unaligned"  & qtype == "Direct"),
                   TRUE,
                   FALSE),
        high =
            ifelse(faction == "Rebels",
                   FALSE,
                   TRUE)
    )

###### combine Assad and Jihad

AJ_Diff <-
    mutate(
        bind_rows(Assad_Diff,
                  Jihad_Diff),
        dv =
            factor(rep(1:2, each = 6),
                   labels = c("Assad",
                              "Jihadis"))
        )

###### plot it

## plot it
plot_both_yeardiff <-
    ggplot(AJ_Diff,
           aes(x    = qtype,
               y    = mean)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Change in Percent Support") +
    scale_y_continuous(limits       = c(-0.40, 0.40),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-0.875, 0.875, .125)) +
    ## misc
    ## scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_grid(dv ~ faction) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_both_yeardiff_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_both_yeardiff),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_both_yeardiff_egg,
       filename = "both-yeardiff-3cat.pdf",
       family   = "Palatino",
       width    = 6.00,
       height   = 2.00)

##### Figure 13: misrepresentation by faction
###### gather data

Dat_Assad <-
    transmute(
        filter(Assad$pool$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct"))
    )

Dat_Jihad <-
    transmute(
        filter(Jihad$pool$main$sub0, # drop "other"
               faction != "other"),
        mean,
        se,
        year = year,
        faction = factor(ifelse(faction == "govt", 1,
                         ifelse(faction == "oppn", 2,
                                3)),
                         labels = c("Government",
                                    "Rebels",
                                    "Unaligned")),
        ## put direct Q above indirect Q in plot
        qtype = factor(qtype,
                       levels = levels(qtype)[2:1],
                       labels = c("List",
                                  "Direct"))
    )

###### process data

## we want the *change* in misrepresentation over surveys:
## (List15 - Direct15) - (List17 - Direct17)
## that's where the c(-1, 1, 1, -1) business comes from below

## ugh, break this up
A_Govt <- subset(Dat_Assad, faction == "Government")
A_Rebs <- subset(Dat_Assad, faction == "Rebels")
A_None <- subset(Dat_Assad, faction == "Unaligned")

J_Govt <- subset(Dat_Jihad, faction == "Government")
J_Rebs <- subset(Dat_Jihad, faction == "Rebels")
J_None <- subset(Dat_Jihad, faction == "Unaligned")

A_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(A_Govt$mean * c(-1, 1, 1, -1)),
                     sum(A_Rebs$mean * c(-1, 1, 1, -1)),
                     sum(A_None$mean * c(-1, 1, 1, -1))),
               se =
                   c(sqrt(sum(A_Govt$se^2)),
                     sqrt(sum(A_Rebs$se^2)),
                     sqrt(sum(A_None$se^2))))

J_Diff <-
    data_frame(faction =
                   factor(1:3,
                          labels = c("Government",
                                     "Rebels",
                                     "Unaligned")),
               mean =
                   c(sum(J_Govt$mean * c(-1, 1, 1, -1)),
                     sum(J_Rebs$mean * c(-1, 1, 1, -1)),
                     sum(J_None$mean * c(-1, 1, 1, -1))),
               se =
                   c(sqrt(sum(J_Govt$se^2)),
                     sqrt(sum(J_Rebs$se^2)),
                     sqrt(sum(J_None$se^2))))

AJ_Diff <-
    mutate(
        bind_rows(A_Diff,
                  J_Diff),
        dv =
            factor(rep(2:1, each = 3),
                   labels = c("Jihadis",
                              "Assad"))
    )

## reorder factors according to means from Assad outcome
A_Diff$faction <-
    reorder(A_Diff$faction, A_Diff$mean)
J_Diff$faction <-
    reorder(J_Diff$faction, A_Diff$mean)

AJ_Diff2 <-
    mutate(
        bind_rows(J_Diff,
                  A_Diff),
        faction =
            factor(faction,
                   levels = rev(levels(faction))),
        dv =
            factor(rep(2:1, each = 3),
                   labels = c("Assad",
                              "Jihadis"))
    )

###### plot it (by dv)

## plot it
plot_misrep_by_dv <-
    ggplot(AJ_Diff2,
           aes(x    = faction,
               y    = mean)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Change in Misrepresentation") +
    scale_y_continuous(limits       = c(-0.15, 0.65),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-0.875, 0.875, .125)) +
    ## misc
    scale_fill_manual(values = c(NA, "red")) +
    coord_flip() +
    facet_wrap(~ dv, ncol = 2) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_misrep_by_dv_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_misrep_by_dv),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_misrep_by_dv_egg,
       filename = "misrep-2dv.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 1.50)

#### Appendix I: Broken Jihadi List
##### Figure 14: broken lists by metric

## gather data
Dat0 <-
    transmute(filter(Jihad$base$main$sub0,
                  qtype == "list",
                  faction %in% c("govt", "none")),
           year,
           metric = factor(metric,
                           levels = levels(metric)[2:1],
                           labels = c("Okay",
                                      "Good")),
           faction = factor(ifelse(faction == "govt", 1, 2),
                            labels = c("Government", "Unaligned")),
           mean,
           se,
           n)

## plot it
plot_jihad_metric_year <-
    ggplot(Dat0,
           aes(x = metric,
               y = mean)) +
    ## vertical line at 0
    geom_hline(yintercept = 0,
               linetype   = "dashed",
               color      = "grey50") +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-.55, .55),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-1.125, 1.125, .25)) +
    ## misc
    coord_flip() +
    facet_grid(year ~ faction) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_jihad_metric_year_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_jihad_metric_year),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_jihad_metric_year_egg,
       filename = "jihadi-metric-1517.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 2.00)

#### Appendix J: Multivariate Results
##### Table 13: list experiment estimates
###### data setup

## for LISTIT, need the list group variable not to include the
## control group's added-values as well
LE <-
    mutate(AC,
           survey17 =
               ifelse(year == "2017", 1, 0),
           assad_list =
               ifelse(ac_group_treat, ac_assad_list, NA),
           caliphate_list =
               ifelse(ac_group_treat, ac_caliphate_list, NA)
           )

###### assad/govt, caliphate/oppn
####### govt

LE0 <- filter(LE, fac_4cat == "govt")

Y <-
    with(LE0,
         cbind(assad_list,
               ac_parl_direct,
               ac_no_parties_direct,
               ac_mil_direct))
X <-
    with(LE0,
         cbind(survey17,
               young,
               male,
               minority,
               lnsrd = log(syr_room_density),
               lnlrd = log(leb_room_density),
               educ_hs,
               interest_bi,
               quran_weekly,
               relpol_bi
               ))
map <- matrix(1, nrow = 4, ncol = ncol(X))
mod <- listit(Y = Y, X = X, map = map)

## gather some estimates
mod_b  <- mod$b.star
mod_se <- mod$se
mod_p  <- pnorm(-abs(mod_b/mod_se)) * 2
mod_n  <- mod$n.treatment
mod_ll <- mod$logL
## spit out something readable
list_b  <- mod$b.star
list_se <- mod$se
list_p  <- (1 - pnorm(abs(list_b / list_se))) * 2
cbind(c("Intercept", colnames(X)),
      format(round(list_b,  2), digits = 2, nsmall = 2),
      format(round(list_se, 2), digits = 2, nsmall = 2),
      format(round(list_p,  3), digits = 3, nsmall = 3))

## gather estimates to shoehorn into texreg
mod_listit_govt_assad    <- mod
mod_listit_govt_assad_b  <- mod$b.star
mod_listit_govt_assad_se <- mod$se
mod_listit_govt_assad_p  <- pnorm(-abs(mod$b.star/mod$se)) * 2
mod_listit_govt_assad_n  <- mod$n.treatment
mod_listit_govt_assad_ll <- mod$logL

## direct Q: logit version
mod_logit_govt_assad <-
    glm(ac_assad_direct
        ~ survey17
        + young
        + male
        + minority
        + log(syr_room_density)
        + log(leb_room_density)
        + educ_hs
        + interest_bi
        + quran_weekly
        + relpol_bi
       ,
        family = binomial,
        data   = LE0)
summary(mod_logit_govt_assad)

####### oppn

LE0 <- filter(LE, fac_4cat == "oppn")

Y <-
    with(LE0,
         cbind(caliphate_list,
               ac_pres_direct,
               ac_multiparty_direct,
               ac_mandate_direct))
X <-
    with(LE0,
         cbind(survey17,
               young,
               male,
               ## minority,                # drop because only 5 people
               lnsrd = log(syr_room_density),
               lnlrd = log(leb_room_density),
               educ_hs,
               interest_bi,
               quran_weekly,
               relpol_bi
               ))
map <- matrix(1, nrow = 4, ncol = ncol(X))
mod <- listit(Y = Y, X = X, map = map)

## gather some estimates
mod_b  <- mod$b.star
mod_se <- mod$se
mod_p  <- pnorm(-abs(mod_b/mod_se)) * 2
mod_n  <- mod$n.treatment
mod_ll <- mod$logL
## spit out something readable
list_b  <- mod$b.star
list_se <- mod$se
list_p  <- (1 - pnorm(abs(list_b / list_se))) * 2
cbind(c("Intercept", colnames(X)),
      format(round(list_b,  2), digits = 2, nsmall = 2),
      format(round(list_se, 2), digits = 2, nsmall = 2),
      format(round(list_p,  3), digits = 3, nsmall = 3))

## gather estimates to shoehorn into texreg
mod_listit_oppn_caliph    <- mod
mod_listit_oppn_caliph_b  <- mod$b.star
mod_listit_oppn_caliph_se <- mod$se
mod_listit_oppn_caliph_p  <- pnorm(-abs(mod$b.star/mod$se)) * 2
mod_listit_oppn_caliph_n  <- mod$n.treatment
mod_listit_oppn_caliph_ll <- mod$logL

## direct Q: logit version
mod_logit_oppn_caliph <-
    glm(ac_caliphate_direct
        ~ survey17
        + young
        + male
        ## + minority
        + log(syr_room_density)
        + log(leb_room_density)
        + educ_hs
        + interest_bi
        + quran_weekly
        + relpol_bi
       ,
        family = binomial,
        data   = LE0)
summary(mod_logit_oppn_caliph)

###### govt + oppn

## models
## NB: will override with listit estimates
## NB: silly to put govt logit estimates in because of
##     absurd SEs due to lack of variation leading to separation
tex_mods <- 
    list(mod_logit_govt_assad,
         mod_logit_oppn_caliph)

## model names
tex_modnames <- 
    c("Government: Assad",
      "Rebels: Jihadis")

## variable names (original order)
tex_vars <- 
    c("(Intercept)",
      "2017",
      "Young",
      "Male",
      "Minority",
      "Crowded in Syria",
      "Crowded in Lebanon",
      "Educated",
      "Interested in Politics",
      "Private Religion",
      "Public Religion")

## override the estimates with listit counterparts
override_b <-
    list(mod_listit_govt_assad_b,
         mod_listit_oppn_caliph_b)
override_se <-
    list(mod_listit_govt_assad_se,
         mod_listit_oppn_caliph_se)
override_p <-
    list(mod_listit_govt_assad_p,
         mod_listit_oppn_caliph_p)
override_n <-
    list(mod_listit_govt_assad_n,
         mod_listit_oppn_caliph_n)
override_ll <-
    list(mod_listit_govt_assad_ll,
         mod_listit_oppn_caliph_ll)

## texreg it to screen
## NB: must replace N and LL with hand-entered versions
texreg(tex_mods,
       custom.model.names = tex_modnames,
       custom.coef.names  = tex_vars,
       override.coef      = override_b,
       override.se        = override_se,
       override.p         = override_p,
       single.row         = TRUE,
       leading.zero       = TRUE,
       booktabs           = TRUE,
       dcolumn            = TRUE,
       stars              = c(.01, .05, .10))
unlist(override_n)
unlist(override_ll)

##### Figure 15: jihadi support among religious and secular rebels
###### gather data

## Assad private religion
le_assad_privrel <-
    sim_me(mod  = mod_listit_govt_assad,
           x    = 10,
           reps = 1000)
le_assad_privrel_pp <-
    calc_pp(sim = le_assad_privrel)
le_assad_privrel_diff <-
    calc_diff(sim = le_assad_privrel)

## Assad public religion
le_assad_pubrel <-
    sim_me(mod  = mod_listit_govt_assad,
           x    = 11,
           reps = 1000)
le_assad_pubrel_pp <-
    calc_pp(sim = le_assad_pubrel)
le_assad_pubrel_diff <-
    calc_diff(sim = le_assad_pubrel)

## Jihadis private religion
le_caliph_privrel <-
    sim_me(mod  = mod_listit_oppn_caliph,
           x    = 9,
           reps = 1000)
le_caliph_privrel_pp <-
    calc_pp(sim = le_caliph_privrel)
le_caliph_privrel_diff <-
    calc_diff(sim = le_caliph_privrel)

## Jihadis public religion
le_caliph_pubrel <-
    sim_me(mod  = mod_listit_oppn_caliph,
           x    = 10,
           reps = 1000)
le_caliph_pubrel_pp <-
    calc_pp(sim = le_caliph_pubrel)
le_caliph_pubrel_diff <-
    calc_diff(sim = le_caliph_pubrel)

## gather the difference estimates
LE_sim <-
    mutate(
        as_data_frame(
            rbind(le_assad_privrel_diff$diff,
                  le_assad_pubrel_diff$diff,
                  le_caliph_privrel_diff$diff,
                  le_caliph_pubrel_diff$diff)
        ),
        group =
            factor(rep(1:2, each = 2),
                   labels = c("Government for Assad",
                              "Rebels for Jihadis")),
        rel =
            factor(rep(1:2, 2),
                   labels = c("Private Religion",
                              "Public Religion"))
    )

## gather the PP estimates, govt
PP_govt <-
    mutate(
        as_data_frame(
            rbind(le_assad_privrel_pp$pp,
                  le_assad_pubrel_pp$pp)
        ),
        val =
            factor(rep(2:1, 2),
                   labels = c("Religious",
                              "Secular")),
        rel =
            factor(rep(1:2, each = 2),
                   labels = c("Private Religion",
                              "Public Religion"))
    )

## gather the PP estimates, oppn
PP_oppn <-
    mutate(
        as_data_frame(
            rbind(le_caliph_privrel_pp$pp,
                  le_caliph_pubrel_pp$pp)
        ),
        val =
            factor(rep(2:1, 2),
                   labels = c("Religious",
                              "Secular")),
        rel =
            factor(rep(1:2, each = 2),
                   labels = c("Private Religion",
                              "Public Religion"))
    )

###### plot oppn

plot_oppn <- 
    ggplot(PP_oppn,
           aes(x = val,
               y = mean)) +
    ## error bars and points
    geom_errorbar(aes(ymin = lo90,
                      ymax = hi90),
                  width = 0.10) +
    geom_errorbar(aes(ymin = lo95,
                      ymax = hi95),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-.05, .80),
                       labels       = seq(-1.0, 1.0, .25) * 100,
                       breaks       = seq(-1.0, 1.0, .25),
                       minor_breaks = seq(-1.125, 1.125, .25)) +
    ## misc
    coord_flip() +
    facet_wrap(~ rel, ncol = 2) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_oppn_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_oppn),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_oppn_egg,
       filename = "le-oppn-rel.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 1.50)

##### Figure 16: jihadi support among rebels

## gather list experiment data
Dat0 <-
    transmute(
        filter(Jihad$pall$main$sub1,
               faction %in% c("fsa", "islamists")),
        mean,
        se,
        faction = factor(ifelse(faction == "fsa", 2, 1),
                         labels = c("Islamists",
                                    "Nationalists")),
        qtype = factor(qtype,
                       labels = c("Direct",
                                  "List"))
    )

## plot it
plot_jihad_oppnall <-
    ggplot(Dat0,
           aes(x     = faction,
               y     = mean)) +
    ## error bars and points
    geom_errorbar(aes(ymin = mean - se * 1.64,
                      ymax = mean + se * 1.64),
                  width = 0.10) +
    geom_errorbar(aes(ymin = mean - se * 1.96,
                      ymax = mean + se * 1.96),
                  width = 0.00) +
    geom_point(size = 1.00) +
    ## percent scales
    labs(x = NULL,
         y = "Percent Support") +
    scale_y_continuous(limits       = c(-0.05, 0.45),
                       labels       = seq(-1.0, 1.0, .20) * 100,
                       breaks       = seq(-1.0, 1.0, .20),
                       minor_breaks = seq(-1.0, 1.0, .10)) +
    ## misc
    coord_flip() +
    facet_wrap(~ qtype, nrow = 1) +
    theme_bw() +
    theme(aspect.ratio = .25) +
    theme(legend.position  = "none",
          strip.background = element_blank(),
          panel.border     = element_rect(color = "black"))

## fiddle with the layout (requires egg package)
plot_jihad_oppnall_egg <- 
    grid.arrange(grobs = lapply(
                             list(plot_jihad_oppnall),
                             set_panel_size,
                             width  = unit(1.5, "in"),
                             height = unit(0.5, "in")
                         ))

## save it
ggsave(plot     = plot_jihad_oppnall_egg,
       filename = "jihadi-poolall-oppn.pdf",
       family   = "Palatino",
       width    = 4.00,
       height   = 1.50)

### (Emacs local variables)

## Local Variables:
## eval: (outline-minor-mode 1)
## eval: (font-lock-mode 1)
## outline-regexp: "^\###+ "
## outline-level: (lambda ()
##                  (cond
##                   ((looking-at "^### ")     1)
##                   ((looking-at "^#### ")    2)
##                   ((looking-at "^##### ")   3)
##                   ((looking-at "^###### ")  4)
##                   ((looking-at "^####### ") 5)
##                   (t 1000)))
## eval: (font-lock-add-keywords
##        nil
##        '(("^### .*" 0 (aref outline-font-lock-faces 0) t)
##          ("^#### .*" 0 (aref outline-font-lock-faces 1) t)
##          ("^##### .*" 0 (aref outline-font-lock-faces 2) t)
##          ("^###### .*" 0 (aref outline-font-lock-faces 3) t)
##          ("^####### .*" 0 (aref outline-font-lock-faces 4) t)))
## eval: (hide-sublevels 1)
## End:
